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SUMMARY 


A tabular method is presented for the determination of the 
spanwise and chordwise bending moments of a rotor blade. The 
method does not require a preliminary smoothing of the mass and 
inertia distributions of the blade, is simple to use, accurate in the 
results, and economical in time. Calculation of the principal 
table, number 4, takes about 3 to 4 days. 

The method is illustrated for the 23-ft. McDonnell blade at an 
advance ratio of u» = 0.35. 


INTRODUCTION 


id IS VITAL to have reliable estimates of the steady 
and alternating stresses to which the blades of a 
helicopter rotor are subjected. Several methods exist 
for the calculation of steady-state stresses, while the 
collocation method," ? recently extended by Yuan* to 
blades with variable mass and moment of inertia dis- 
tributions, is widely used for determining the alternat- 
ing stresses. 

The tabular method presented in this paper is an 
adaptation of the method introduced by Myklestad 
in 1944.4 Except for the one serious drawback that 
the unknowns of the problem are determined as small 
differences of extremely large numbers, so that the ex- 
isting ten-row calculating machines are barely adequate 
to handle the majority of problems, and are inadequate 
to handle some others, the method has numerous dis- 
tinguishing features: 

(1) It is a numerical method. As such it is appli- 
cable to blades with arbitrary mass, chord, and moment 
of inertia distributions. (When the collocation method 
is applied to the analysis of a blade, the distributions 
must first be “properly” smoothed out, otherwise the 
convergence of the method will be exceedingly slow and 
the results may be considerably in error.) The blade 
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* The author wishes to express his indebtedness to his col- 
leagues, Elizabeth J. Spitzer, Kathryn Meyer, and Frances 
Schniitz for carrying out the calculations and for checking the 
derivations. 
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is subdivided into small sections and the forces are 
assumed to act at the station points‘of the sections. 
These need not be the midpoints of the sections. 

(2) The amount of work involved is proportional to 
the number of stations used. (In the collocation 
method the amount of work increases exponentially 
with the number of stations.) Sections of unequal 
lengths are handled with no more effort than sections 
of equal lengths.f 

(3) The determination of the afternating stresses is 
similar to the determination of the steady-state stresses. 
No new principles are involved; only additional tabu- 
lation. 

(4) The method is not one of successive approxima- 
tions. The unknowns of the problem are the coef- 
ficients Bp, 2, . . . . of the blade tip slope 


ip = = + Bcosy + & sing + 


(39a) 
and the coefficients RY:, RY¥2, .... of the 
blade tip deflection 
Nip = RY = RYO + RVicose + Rosin g + 

RY3;cos2¢+.... (39b) 


Working station by station from tip toward root, the 
station bending moments, slopes, deflections, inertia 
loads, air loads, centrifugal loads of the various har- 
monics are computed in terms of the tip @, and Y;. Im- 
position of the boundary conditions at the hinge per- 
mits determination of the unknown ©, and Y,. 

(5) Problems with homogeneous boundary conditions 
(such as vanishing deflections and bending moments 
at the origin when the blades are hinged), or with more 
complicated boundary conditions (which arise, fur in- 
stance, when two rigidly connected blades are mounte«d 


t No use of this fact was made in the analysis of the 23-ft. 
McDonnell blade. However, more recent calculations take ad- 


vantage this feature. 
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in seesaw fashion, or when there is a blade damper, as 
in lag motion, etc.) are handled with equal ease. 

(6) The effects of the various load components are 
kept separate and are readily interpreted. The steady 
and alternating hinge forces are obtained without sepa- 
rate calculations. 

(7) The method permits elimination of the assump- 
tion, vital in all previous studies, that the air forces re- 
main unchanged in the course of the flexing of the blade. 

(8) The method is in principle exceedingly simple 
(the bending vibrations counterpart of Holzer’s well- 
known method for torsional vibration calculations), and 
skill in performing the calculations can be readily ac- 
quired by computers: unfamiliar with engineering ter- 
minology. 

The development of the tabular method was under- 
taken in the spring of 1946,.under sponsorship of the 
Bureau of Aeronautics, U.S. Navy, when it was realized 
that the slow convergence of the collocation method 
(as the number of stations was gradually increased 
from five to eight) cast doubt upon the correctness of 
the method of approach.” This, coupled with such well- 
known drawbacks of the collocation method as the 
necessity of using extremely smooth distribution func- 
tions that blur out the effect of any sharp changes in- 
troduced in the blade design, made the search for a more 
satisfactory method imperative.* 

In addition to th@novel method of bending moment 


calculations, the present paper contains a number of. 


minor results which, to the author’s knowledge, were 
not heretofore available in the literature. These are: 
an expression for the air force in forward flight, Eq. 
(37), which takes into account the elastic deflection of 
the blade; an expression for the air force, Eqs. (69), 
(70), which takes into account the lag oscillations; an 
equation, Eq. (80), with appropriate boundary condi- 
tions, Eq. (72), for the elastic deflection of the blade 
in the plane of rotation; and a formulation of the 
boundary conditions, Eqs. (67), for teetered blades. 

A few weeks ago, after completion of most of the 
present work but before the final write-up, eight excel- 
lent reports by W. C. Johnson, R. Mayne, and others 
of the Goodyear Aircraft Corporation’ came to the at- 
tention of the author. The Goodyear reports deal 
with the same problem dealt with in this paper; 
namely, the establishment of a satisfactory tabular 
schedule to replace the collocation method. The tabu- 
lar method developed by Goodyear and the one de- 
veloped herein are identical in essence, but unlike in de- 
tails. The Goodyear group also developed a quasi- 
static method for the solution of the dynamical prob- 
lem. This is an excellent short-cut procedure which 
the author also attempted to obtain, but without 
success. 


* For instance, one of the puzzling questions—how does a con- 
centrated mass placed into the tip of the blade affect the bending 
moment—cannot be treated by the collocation method except 
in a slipshod fashion. 
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As an illustration of the method, the McDonnelj 
fabric-covered 23-ft. blade (NACA 0019 section at 
r = 33.5 in. tapered to 0013 at tip) is analyzed for 
bending moment in flapping and in lag (with a lag 
damper attachment of 400 Ib. ft. per rad. per sec.) at 
u = 0.35, and the procedure for stress analysis is also 
given. The spanwise bending moments are deter- 
mined by a 13-station tabular calculation, by an g. 


‘station collocation calculation, and by a 13-station 


Johnson-Mayne quasi-static calculation (adapted to 
the procedures of this paper). The results are in ex. 
cellent agreement. 


NOTATION 
R = rotor radius 
r = distance of blade element from axis of rotation 


(flapping hinge distance from center of ro- 
tation considered negligible) 


x = 7/R 

m'(x) = blade mass per unit length at station x 

ma(x) = 2(1 — see Eq. (80) 

c(x) = chord length 

I(x) = cross-section moment of inertia of blade spar 

E = modulus of elasticity 

w = rotor frequency 

t = time 

v. = wt = blade azimuth (counted from rearmost 
position of blade) 

b = number of blades per rotor 

r = inflow ratio 

m = advance ratio 

WwW = weight of helicopter 

n = load factor, see Eq. (19a) 

P = shaft horsepower of helicopter (includes factor 
of efficiency) 

m = power factor, see Eq. (19b) 

p == air density 

v =: blade pitch angle 

xt+de = angle of attack 

cL = a(x + 8) = lift coefficient . 

cp = Ko + «2(x + 8)? = coefficient of profile drag 

U = resultant air velocity at blade section 

Up = perpendicular component of air velocity 

Ur = tangential component of air velocity 

Ri = length of lag hinge link 

u = (r — Rx)/(R — Rt) 

Ms = Ser m'ar.= total blade mass 

Rec. = Ms  m'rdr = c.g. distance of Ms from 

center of rotation 

Ir = Sor m'rtdr = moment of inertia of Mg with 
respect to flapping hinge, see Eq. (+4) 

Ms = JSgim'dr = outboaril blade mass 

= Mi Spam'(r — Ri)dr = c.g. distance of Ms 
from lag hinge 

It = — Rz)*dr = moment of inertia of My 
with respect to lag hinge, see Eq. (76) 

€ = MRiRe.g./I, = lag hinge eccentricity para- 
meter, see Eq. (76) 

D = damping constant of lag damper 

K = spring constant of centering spring (also, 
stress concentration factor) 

d. = length of blade section k 

le, ett = distance of station k from station k + | 

Nk+t = distance of outboard end of section & from sta- 


tion k + 1 


: 
i 


STRESS ANALYSIS 


g,h, u,v = deflection and slope coefficients, defined by 
Egs. (9) 

Tq ' = distribution functions involving mass and 
chord, see Eqs. (1g) and (16) 

C(x) = centrifugal tension at station x 


6S, 6J, 6G, 6Cr, 6T, 6C, = increments of various forces, 
cf. Eqs. (1) and (69) 


bL st. = static lift force expression, see Eq. (18) 

bLex. = exact lift force expression, see Eq. (37) 

5Lrig. = lift force acting on rigid blade, see Eq. (41) 

5L gs. = Johnson-Mayne’s quasi-static lift force, see 
Eq. (58) : 

6D = '/:pccpU*dr = drag force, see Eq. (69e) 

r = average thrust per blade 

= Sr = average torque per blade, see Eq. 
(40b) 

= — = torque with respect to 
the lag hinge, see Eq. (75) 

cr = 6T/xpw*R‘ = thrust coefficient, see Eq. (30) 

y(r) = + n(r) + sin + cos 
2y + .... = deflection of blade (flapping 
or lag) 

{y}e = deflection at station k 

T(r) = 

= dy/dr 

¢'(r) = dy/dr 

Y = Juip/R 

= ¢tip 

S(x) = So(x) + Si(x) cos ¥ + S2(x) sin y + S;(x) cos 
2y + .... = shear force [flapping or lag; 
strictly speaking — S is the shear force, see 
Eq. (2a) ] 

M(x) = bending moment (flapping or lag) 

Mr, Mi = spanwise (flapping) and chordwise (lag) 
bending moments (when such distinction is 

necessary) 

B = +a, cosy + bh siny + a:cos2y+.... 
= flapping angle of rigid blade 

= Ayo + Aicosy + B,siny + Ascos 2y +.... 
= lag angle of rigid blade 

a, b, c, d (with subscripts M, S, etc.) = coefficients in Eqs. (11) 
and (33) 

Tr, Um, etc. = moments of the distribution functions r, g, see 
Egs. (17) 

1 = angles in Fig. 6 


cos ¥ — components of tip slope for various 


= 
: conditions, see Eqs. (53), (54), (55) 


S* = cos ¥ — component of auxiliary shear force, 
Eq. (55) 

CG = coefficients in the collocation expansion of T, 
see Eq. (83) 

P. = polynomials in the collocation method, see 
Eqs. (82) 


a = 0, a1, a2, as, ay = coefficients defined by Eqs. (82) 


QT, Ty = amplification factors, see Eq. (59) 

» = natural bending frequency of blade 

sin rx — H» = deflection curve of hinged blade in first natural 
mode 

E., E.,E, = kinetic, centrifugal, and strain energies, see 
Egs. (61) 

fe = bending stress 

fi = tensile stress 

Zz = section modulus 

A = area of cross section 


OF ROTOR BLADES 


Fi, Fi, Fu, Fo, fmoz., fmin., fm, fe = stresses defined in the sec- 
tion on stress analysis 


Ro, R: = ratios defined by Eqs. (85) 
MSuit. = ultimate margin of safety 
M Stat. = fatigue margin of safety 


OUTLINE OF THE TABULAR PROCEDURE FOR HINGED 
BLADES IN FLAPPING 


The basic assumptions of the method are: the blade 
is torsionally rigid; lag motion does not affect the flap- 
ping motion (but conversely it does); inflow velocity 
and lift coefficient are constant through the rotor disc; 
the thrust component of air drag and the effect of radial 
air velocity are negligible; controls are in neutral posi- 
tion; deflections and angular changes are small, their 
higher powers can be neglected. 

Fig. la shows the flexed blade, subdivided into nine 
whole sections, d;, dz, . . . ., dy, plus two half-sections 
d,,, and dg,. Two zero sections dy and dy are added 
for convenience. In the tabular method one replaces the 
given stiffness distribution, mass distribution, and plan 
form of the blade by step-function approximations EI(x), 
m'(x), and c(x) which have the constant values { E/},, 
m,' and ¢,, respectively, over each section k. Also, 
the forces that act on the sections are considered as con- 
centrated into the station points of the sections.* 


* means: the slope gat station k; { ¢} (or{¢} means: 
¢ at the right end (or left end) of section dy. [The symbol gx, 
without braces, is reserved for another purpose, see Eq. 
(36).] The notation for the functons ¢’ (curvature), y (de- 
flection), M (bending moment), and S, 5L, etc., (forces) is similar. 
On the other hand, the prescribed constant m'd at station k may 
be written either as mz'd, or as{m‘d} x, according to convenience. 
In place of station & it will be convenient to refer occasionally to 
station x or perhaps xy. The notation then is, or{ 


Fic. 1. Diagram (a), flexed blade. Diagram (b), equilibrium 
of blade section S are positive upward, ¢’ and 
positive for compression of upper fiber. 


| 
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The blade element {m'd},, drawn enlarged in Fig. 1b 
for k = 5, has to remain in equilibrium under the action 
of the applied forces 

éL = lift force, see Eqs. (18), (37), (41), (58) (1a) 
—éG = dead weight = —{m’dg}, (1b) 
= inertia force = —{m'dy}, (1c) 

shear component of centrifugal force = 
= 


éT = shear resultant of tension force = 
{oe} er (le) 


= 


(1f) 
where 
tip 
t, = Rw*{m'd},, C, = centrifugal tension = Lire}. 
(1g) 
and the reacting shear, —5S, thus: 
0 = + — + BL — + (2a) 


Note that the drag and the radial components of the 
air force, of 6J, etc., have small effect, and can be con- 
sidered as negligible. Note also that the formula relat- 
ing slope and curvature 
{e}er — 

is approximate, but adequate for short sections. 

A moment equilibrium of cement {m'd}, must also 
exist: 


(3) 


0 = + (2b) 


Here /, .+1 denotes the distance of station k from sta- 
tion k + 1. 

Hence, if the bending moment and shear are known 
for station k + 1, they can be determined for station k 


from £ 
= {Mp t+ {5M}, = {Mp the, (4) 


{She (Shee 8G + — + (5) 


But at the tip of the blade the bending moment and 
shear force vanish, 


Egs. (4) and (5) thus furnish an algorithm for deter- 
mining M and S step by step, from tip to center of rota- 
tion. 

However, at any station & not only {M}, and {S},, 
but also{ ¢’ {e},and{y}, must be determined. This 
is because {5S}, depends explicitly on curvature, slope, 
and deflection at station k. This renders the present 
calculations more cumbersome than the static calcula- 
tions met in most ordinary analysis. 

{o’}, is readily determined. It is 
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= {M/EI}, (7a) 


{ye}, and {y}, are obtained in a somewhat more compli- 
cated manner from 


fy}e = + Loy}, 


(7b) 
(7c) 
where 


(Sb) 


Here 


Elvi EI, 


is the slope change at station & due to a unit moment 
acting at station k + 1; 


(9a) 


Ur, = 


2EI, 
is the slope change at & due to a unit shear force at k + 


1; and also the deflection change at & due to a unit | 
moment at k + 1; and, finally, 


(9b) 


Ux = = 


Le ct? — 
3EI, 


is the deflection change at station & due to a unit shear 
force acting at station k + 1. (When station k + 1 
is at the’center of section k + 1, then +1 = '/sde+1.) 
Eqs. (9a, b, c) are the well-known slope and deflection 
formulas for a cantilever beam loaded at the tip. 

The difficulties that attend the present calculations arise 
from the fact that, unlike the tip shear and the tip bending 
moment, the tip slope 


= (9c) 


Vip (10a) 


and the tip deflection 

Nip = RY (10b) 
do not vanish. This necessitates carrying along the quan- 
tities @ and Y as unknowns in ail computations. For 
instance, at station x the quantities M, y, y, and S are 
obtained in the form 


Fic. 2. Bending moment of 23-ft. blade at u = 0. 1g load, full 
power condition. 


— | 
Muy = 0, Stip = 0 | 
| 
KONE 
| | 


STRESS ANALYSIS OF ROTOR BLADES 


#0 
— 
——| 

Ti/d © airlead/f+ 
ld = deadusigit/ ft 

2 4 40 


Fic. 3. Loads per foot run acting on blade at » = 0. 


M(x) = ay(x) + by (x)® + Cu(x) Y 
=a,(x) +b,(x)® + o,(x)¥ 
y(x) = a,(x) + b,(x)® + o(x)¥ 
S(x) = a,(x) + b,(x)® + c,(x) 


where the a’s, b’s, and c’s are numerical coefficients, de- 
termined by the step-by-step procedure, but the $ and 
Y are unknown parameters. 

The process is continued until the origin, x = 0, is 
reached. There the two equations (boundary condi- 
tions) of zero moment and zero deflection 


* = 0)y(0) = 0 


permit numerical determination of the two parameters 
@and Y. Resubstitution of these values of 6 and Y 
into the expressions M(x), S(x), y(x), and g(x) furnishes 
the bending moment and shear load distributions, as 
well as the deflection and slope of the blade. 

It will be convenient to refer to the calculation of 
M, ¢, 9, S at station k by Eqs. (4), (7b), (7c), (1), and 
(5) as the long schedule, because it involves the deter- 
mination of the load increments —éC,, 57, etc., prior 
to the calculation of the shear. However, determina- 
tion of { S}, can be carried out more conveniently, with- 
out the intermediate step of computing loads, by the 
formula 


(11) 


(12a, b) 


= {S, + So + Ss + Sole (13) 
where 
[{8L}, is given by Eqs. (18), (37), (41), (58)] 
= {So}ati + {-8G}, (14b) 
[{ 8G}, is given by Eq. (1b)] 
= + (14¢) 
[{8J}. is given by Eq. (1c)] 
{Sole = — {Cor} (14d) 


are the air force, dead weight, inertia, and centrifugal 
shears. 

The calculation of S,, Sg, S, does not differ in essen- 
tials from the calculation of Eqs. (la), (1b), (Ic). In- 
stead of tabulating an increment, one tabulates at once 
the increment plus the entry of the preceding station. 


=< 
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Fic. 4. Shear forces acting on blade atu = 0. Z4S5L = air- 
force shear, £ — 6G = deadweight shear, 2{—éCr + 67} = 
centrifugal shear, 25S = total shear. 


However, the step (14d) constitutes a genuine improve- 
ment in the method, inasmuch as it replaces three steps: 
calculation of —5C,, of 5T, and their summation from 
station & to tip, by a single step, calculation of C;,{ ¢},. 
In general, one therefore will use the short schedule, 
Egs. (4), (7b), (7c), (14), in bending moment analysis, 
unless the appraisal of the component loads —5C, and 
5T is the specific purpose of the calculations (as in Figs. 
3 and 4); in this case the long schedule calculations 
must be used. 


BENDING MOMENT OF THE MCDONNELL 23-FrT. BLADE 
AT p= 0 


The fundamental information in regard to the Mc- 
Donnell 23-ft. blade, as needed in this analysis, is con- 
tained in the values 

b = 3 (two rotors) } 
W = 11,000 Ibs. 
R = 23 ft. 
w = 23.35 rad. per sec. 
P =0.85 X 900 = 765 hp. , 
a@ =5.75 
x = 0.0075 
& = 0.8 
= 29 X 10° Ibs. per sq.in. 
p = 0.002378 slug per cu.ft. | 


in the mass distribution m’(x) and moment of inertia 
distribution /(x) tabulated as step functions in Table 1, 
and in the blade chord function c(x) which is given as 0 
from r = 0 to 33.5 in., as 1.25 ft. from 33.5 to 262.2 


(15) 
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OF 23 Fr Biape 
in, 
2176 27% |4010 o 
27¢ 
269-1 |.995\ | 1360.59| | 126.827 0900 
262.2 4725|-8/0 3426943 
1-9 12.3) | 2944.80 | 3777/47) 58 
2346 -85294/# OFF 4327729/F¢0 
220.8 29949 2869.52 
207 | 6634003572 10954180797 
493.2 |-7 a 28/7. 02 
165.6 4 09632 2778.09 |97//.5 31 \2256 
424.2 « | | | 
110.6 1020 ” 294F/.9/ 
96.6 a | «37359386 
82.8 |.3 1038 2993.03| /3228.559\ 3435, 
SF.2 |.2 | .s807 ” 521. 804270.9/9 
“ -220€66/6773 |. ZB525860 33/5 
27.6 4 12-3) | 1622.28) /3922-4 | 59 
6-9 .8060 
-0328/8 07372 | | -0036/6025208 
| o o o 5953.72! |.8700 


in., as 1.0 ft. from 262.2 to 271 in., and as 0 from 271 
to 276 in. (simulating tip loss). 

The column d of Table 1 shows the width of the sec- 
tions, the column / the distance between station points. 
Column 7; is Rw*{m'd} x column C, (the centrifugal 
tension) is 24” {xr},. The air-force column 


de = cad}, (16) 


takes into account tip loss and blade root location at 
33.5 in. (Aerodynamically section d(0.975) extends 
from x = 0.95 to x = 0.98, and therefore d(0.975) 
0.3 X 2.3 is used for q(0.975); similarly, d(0.1) 
0.2884 X 2.3 is used for g(0.1).) 

On the basis of Eqs. (9) one calculates the slope and 
deflection coefficients v, u, h, g to ten digits. They 
constitute the last three columns of Table 1. 

From the 7 and gq distributions one obtains the mo- 
ments 


ty = C(O) = D?xy7, = 15,953.721 Ibs. 
Ty = = Tpw?/R = 9,517.055 Ibs. 
Qo = Uae = 48,335.990 Ibs. 

= = 26,433.019 lbs. 

= = 17,409.991 Ibs. 

= = 12,737.037 Ibs. 


The commonly accepted expression for the air force 
(acting on section d,) in vertical flight is 


= qe x + 


(17) 


(18) 


For untwisted and torsionally rigid blades, as they 
are usually considered in bending moment calculations, 
one has #(x, ¥) = constant. Usually, also the inflow 
ratio A(x, w) is considered as constant, although in 


reality it depends both on the distance 7 of the element 
from the flapping hinge and on the azimuthal location 
y of the blade in the plane of rotation (due to interfer- 
ence effects from the fuselage, pylon arms, and perhaps 
a second rotor). A variation of \ and # with x offers 
no difficulties in the calculations, while a variation 
with y can be treated along the lines developed for 
p> 0. 

It is customary to determine \ and 3 (when constant) 
from the two conditions (note that 6 = 3 and there 
are two rotors) 


T = nW/6 (19a) 
= 550mP/6 (19b) 
where 
n = load factor, m = power factor (20) 
tip 
T= = gv + qn? (21a) 


is the average thrust (over a revolution) per blade; and 


Q =F 00 = + + anX 
(2x2 — + gi(xe — a)A*] (21b) 


is the average torque per blade. 
For the 23-ft. blade one thus finds from Eqs. (19) at 
normal conditions, n = m = 1, that 


T = 1,833.333 Ibs. 
Q = 3,003.212 ft.lbs. 
and, therefore, by Eqs. (21) 


(22a) 
(22b) 
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TABLE 2. Long CALCULATIONS 


For 28 FT BLADE. M20, 


Litem | coffscient of 1 |Coeffictent of 
mM ° 
167.663 167.66 
-46 -2.968 
“M2L.827 200.52 
dr 
|-35.83 
284.098875 |\-6/.8/ 
575.875 575.88 
-46 -7.562 -7 
2730F 7072 147789 
or FOF 724 BAIS |-23.62 
BEVSTE | 3/.57 
2279. 579797 10.81 
$30. 1/30 $30.43 
-7.3¢8 -7.37 
dr | #72 
FASC CEI \-2408S 397.20 | 327063 
| 675-8525 3245 |./237/ 
-5.086 
SG -35.749 
SG (6727/5002 |-/72-23 
dr 4296/05 |-728279/.2F2 | 
2587355.838 |-/a9-56 | 
104.50 
AS SSS TI2 H3BC-67/093 «/2078 
-IG -297.8#6 
SEIF3-SEUF \-336300-7062 |+35:/0 
dr 9/6 210.0697 | 1.0 
Ss e356 2279-63F |-200/88 93.0F |-/83-F8| 
Fro 
IG 
o 


A = —0.04937849 


(23a) 
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Skipping the items 
Mtip = 0, = ® (25) 


one tabulates in Table 2, by the “long schedule,” in suc- 
cession 


M(0.975), (0.975), 6L(0.975), —#éG(0.975), 


—6Cp(0.975), 67 (0.975), S(0.975), M(0.9), .... (26) 
and so on, until the moment at the origin, 
M(0) = 10,797,946.70 — 60,677,944.816 (27) 


is reached. (For illustration of the method it was 
deemed sufficient to reproduce in Tables 2, 3, and 4a 
the entries for stations 0.975, 0.9, 0.8, 0.1, 0.025, 0, and 
omit those for stations 0.7, 0.6, 0.5, 0.4, 0.3, 0.2.) But, 
because of the hinge at the origin, M(0) = 0. This 
gives 


Ytip = & = 0.1779550500 (23) 


and the numerical values listed in the Results column 
now follow atonce. Sand its components are obtained 
in Ibs., M in Ib. ft., and ¢ as dimensionless number. 
M(x) is plotted as curve A of Fig. 2; the loads per ft. 
run, 6L/d, —iG/d, —5C,y/d, 5T/d, and 6S/d are plotted 
in Fig. 3; the shear forces Dt? 6L, =” — 4G, ri? — 
and 5S in Fig. 4.* 

* Figs. 1, 6, and 11 of this paper illustrate general principles 
and are applicable to any blade. The other figures refer to the 


McDonnell 23-ft. blade and are drawn for n = m = 1, except Fig. 
5 in which the load factor is varied from '/; to 3. 


3. Sworr For 22 FT. 


3 = 0.18027326 (23b) 


As a result, the air load has the values 6Z entered in 
column 1 of Table 2. The dead-weight loads —éG, 
Eq. (1b), are also readily computed, and are likewise 
entered in column 1 of Table 2. Since at u» = O the 
air force is constant in time, there are no reacting inertia 
forces, 


6J + 0 (at » = 0, A = constant) (24) 


and the tip deflection RY does not explicitly enter the cal- 
culations. Eqs. (11) appear correspondingly in simpler 
form, without the item y(x), and without the coef- 
ficients Cy, Cp, Ch 
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_ Table 3 shows the bending moment calculations by the 

“short schedule,” Eqs. (4), (7), (14), leaving the inflow 
ratio \ and the blade angle # of the air-force expression 
6L unspecified. This is an advantage when a number 
of load and power conditions are to be investigated. 
Carrying the additional columns } and #, one finds the 
hinge moment in the form 


M(0) = 29,543,801.55\ + 26,000,417.860 — 
51,739.83097 — 17,717,936.75@ (28) 
For 
= —0.04937849 (29a) 
8 = 0.18027326 (29b) 


‘(load factor = power factor = 1), the condition M(0) 
= 0 yields 


® = 0.1792879156 (29c) 


and a “Results” column can now be computed as in 
Table 2. The bending moment based on Eq. (29) is 
plotted as curve B of Fig. 2. The slight deviation of 
the curve A from B arises from the finiteness of the 
sections d,. 

Other load and power conditions are treated similarly. 
The m = 1; n = '/2, 1, 2, 3 curves are plotted in Fig. 
5a, and the corresponding #, X, ® in Fig. 5b. It is 
seen from Figs. 5 that when the load factor is less than 
1, then the helicopter is rising, 

= —0.04912 (30) 
and ‘the maximum bending moment remains roughly 
unchanged in magnitude but moves toward the blade 


tip. When the load factor is greater than 1, then the 
helicopter is sinking, the maximum bending moment 


A> Nnovering 


Fic. 5a. Variation of bending moment with ioad factor for full 
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Fic. 5b. Variation of inflow ratio, blade angle, tip slope with 
load factor. 


remains at about 20 per cent span but increases in 
magnitude with the load factor. 
The writer also calculated from Table 3 the bending 


moment curves m = 1, nm = 0.75, 1.0, 1.25. These . 


curves are similar to the m = m = 1 curve, with peaks 
at x ~ 0.2, of relative heights 93 per cent, 100 per cent, 
107 per cent; they do not exhibit any novel features 
and are not reproduced here. 

The reader by now has probably noted that ten digits 
(the full keyboard of a calculating machine) were car- 
ried in all computations. « Is this really necessary, or is 
something less also sufficient? 

The import of this question can be seen from the 
following. With the values given in Eqs. (29a), (29b), 
and (29c) one obtains > 


M(0.025) = 84.7 (31a) 


while using Eqs. (29a) and (29b) with a five-place ac- 
curate tip slope 


® = 0.1792900000 (32) 
the incorrect 
M(0.025) = 55.1 (31b) 


results. If the step function m’(x), I(x), and c(x) dis- 
tributions of Table 1 are accepted as exact, and 


M(x) = adu(x) + + + dy(x)3 (33) 


is to be determined to 1 Ib. ft. accuracy, then a,(0) 
must be determined accurately to 1 Ib. ft., bj,(0), 
¢u(0), dy(0) to 10 Ib. ft. each (since #, \, are each 
about 0.1 in magnitude). The argument is the same 
for S(x) if this is to be determined to 1-lb. accuracy. 
The author repeated the calculation of column & of 
Table 3, determining the © coefficients, (x), consist- 
ently to two decimal places for M, S, and its compo- 
nents, and to one more digit for g than for M. It was 
found that by(0) = —17,717,944.25 as compared with 


{aL 


and 


the e 
j 


> with 
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the. correct value —17,717,936.7. This indicates that 
an error at the tip station spreads about four places to 
the left by the time the origin x = 0 is reached. Evi- 
dently, therefore, some digits can be omitted in the 
calculation of the outboard stations. Nevertheless, 


it is best practice not to omit any but to carry consist- 
ently ten digits in all calculations.* 
* On the other hand, primary data, such as rz, ge, Jz, need not 


be given accurately at all. But quantities derived from these, 
such as u, v, etc., must again have full keyboard accuracy. 


THE AIR-FORCE EXPRESSION AT p > 0 


When the helicopter advances at the rate u, then the air force® ” 


{aL}. '/sp{cc,U"d} = UpUr + 8Uz*};/w?R? = 


+ cos + sin + SLs cos + bly sin Wi, (34) 


becomes a harmonic variable of y, because the perpendicular component of the air velocity 


Up = [A — ¢(x)u cos — y(x) 


and the tangential component of the air velocity 


Ur = [x + usin YJoR 
Correspondingly ¢ and y also become harmonic variables. Writing all expressions con- 


are harmonic variables. 


(35a) 


(35b) 


sistently to second harmonics, and neglecting higher components, one obtains, with 


o(x, ¥) = go(x) + ¢i(x) cos + sin + cos 2b + g(x) sin 2 
y(x, ¥)/R = T(x, = To(x) + Ti(x) cos + sin + Ts(x) cos + Ty(x) sin 2y 


the exact lift force expression 


(36a) 
(36b) 


= [x + P(x? + + — — Qa] + 
[—uxgo — xT2 — + — cosy + + Quxd + xT, — + — sin + 
— — — cos + — — + 2xTs] sin 2}, (37) 


and the inertia, dead-weight, centrifugal, and flexure forces 


{oJ}, = {Ti cosy + T sin y + 47; cos 2Y + 47, sin 2p}, (38a) 

— {8G}, = —{m'dg}, (38b) 

= + gr cosy + sin + cos WY + sin 2}, (38e) 
(38d) 


= {Cd}a loo’ + + sin + gs’ cos 2p + yy’ sin 2p}, 


The bending moment calculations can be carried out according to the long schedule or the short schedule, as in 
the preceding section, but it must be borne in mind that the tip slope 


=6 = 9% + %, cosy + & sin y + 4; cos 2y + sin 2y 


and the tip deflection 


Yip = RY = RY¥o + RY: cosy + R¥2 sin y + RY; cos 2p + RY, sin 2 


(39a) 


(39b) 


now involve five components each. Thus, in place of the two columns, 1 and 4 of the static calculations, one now 


has to deal with ten columns: coefficients of lo, Bo, 11, By & Be, Y; & Vo, 1s, 13, Bs & Ms, V3 & Vs, la. 


Note 


that Y, does not explicitly enter any of the formulas and therefore its column can be omitted, and that, furthermore, 
the coefficient of &, is the same as the coefficient of ,, the coefficient of Y, the same as the coefficient of V3, etc. 
In the air-force expression, Eq. (37), \ and # are as yet unspecified quantities. They again can be determined 


from Eqs. (19a, b), where riow 


T= = dx + (x? + + — — 


(40a) 
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ti 
Q= = (R/a)d {xg} (mo + (x? + + — + — — + 


— — Apes + 1/281? + 1/282? + 283? + + + — — + Tai) + 
+ + + + + (40b) 
(8Dp is the constant part of the lag component 4D of the air force, cf. Eq. (69e).) But these two equations cannot 
be used for calculating \ and 3 without prior knowledge of g(x) and y(x). It is therefore best to perform the bend- 


ing moment analysis for assumed values of \ and # (for instance, by assigning \ and # their “rigid blade” values, 
see further below), and then, after the g(x) and y(x) functions have been established, verify the agreement with or 


departure from Eqs. (19a) and (19b). 
Because of the extreme tediousness of calculating with 5L.,., it is usually preferable to use the rigid blade lift 


force approximation 
— (x? + + cos + [urd + + (x? 
— 
which is obtained from 6L,,, by writing 
¢=T=B8 =a +a, cosy + b sin y + a2 cos + sin (42) 
in Eq. (37). This approximation is adequate in most cases.* 4 is called the flapping angle of the rigid blade; 
4, 1, by, ...., the flapping coefficients. (Watch the positive signs in the expansion of 8; this is contrary to the 
conventional notation. ) Correspondingly T and Q must be written in the form 


T = gd + (gu + — (43a) 
Q = (R/a)[(Ko + (Quix + + — — + 
— a) — wrar) + + + + + gru(aobs + — + 
+ + + + + '/20002)}] (43h) 
d and # (and the flapping coefficients a, ...., 62) can now be predetermined from Eqs. (19), (43), and the five 
equations into which the condition of moment equilibrium with respect to the flapping hinge, 


— + sin 


IpB + Ipw*B = —MagRoc. + pce, U*rdr (44) 
decomposes: 
cos 0: T1140 = + + — — 
cos yp: O = — + + 
sin y: O= + 2ugu? + — + (45) 
cos 2: —3rya2 = — — 
sin —3ryb2 = —'/ou*qido — + 


* The Goodyear group determined experimentally the bending moment of a hinged blade, simulating both éZrig. and 5Lex. condi- 
tions, and found good agreement between the two results. 


— 2x*%b.] cos + [—1/2u%a9 — + sin 2}, (41). 


BENDING MOMENT OF THE MCDONNELL BLADE AT 
pw = 0.35 


For the McDonnell 23-ft. blade the “rigid blade”’ 
equations, Eqs. (19), (43), and (45), lead to 


& = 0.21464752, A = —0.09572523, a = 0.14385868, 
a = —0.14659538, b} = —0.06780956, a = 
—0.01377962, = 0.00598001 (46) 


and the air force assumes the form 


{ }e = Ge { [0.01296402 — 0.09572523x + 
0.21464752x?] + [0.00207667 — 0.05276197x + 
0.06780956x?] cos y + [—0.02901435 + 0.15129977x — 
0.14659538x*] sin y + [—0.01314716 + 0.05130838x — 
0.01196002x*] cos 2¥ + [—0.00881134 + 0.02373335x 
0.02755924x?] sin (47) 


A summation, S, = >”, yields the air-force shear at 
station k which is entered, along with the dead-weight 
shear Sg, as the Df? 8G} item of colums lo, hy, 
le, 13, 1, of Table 4a. (The following checks are avail- 
able: = 7, = = 0, = 
= The remainder of Table 
4a is completed in accordance with the short schedule. 
(For the time being disregard the column 5;* and 5;*.) 

The boundary conditions 


= 0, = 0 (6) 


reduce to the nine equations 


M,(0) = 3,143,209.682 — 17,717,936.75%) = 0 (48a) 
M,(0) = 335,653.5109 — 24,223,138.93(4, — Vi) = 0 
(48b) 
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1 
TABLE 46, Resu7Ts of ms BENDING MOMENT CarcuLATIONS ATM.3S, 
— mi 
Tren | cosO | cosy |swy |cosay|sway | Coso | Cosy] |CoS2¥ sway hat 
inboard Sebati 
fe | taboard selution 
mM 77. 1064-29 | -1/0.68| 108./2| ~2.8/ tha 
SL-IG | 105.3) | 20.5¢|-28.37| 34.71 |-1.15 | 1099.72 | 12.78 | 32.05 | 507,20 
Fo |-36.93 | 26.79 | | 106-48 |- 558.98) 267.77 Fig 
S999 | 18.86|-3.57 | -6.77 | -CY -1693.25 | |-/2/.53 | 6¢.84 - 74-93 (wh 
| ¢.52 | 1-780 | 3.87 Ss 6-47 |-2.586 | /70/ | 13-57 | 3.34 nor 
slog 
mM | 7.80 |-/659|-/3.45 | 6.68 74 4297.17 | | - 139.32 | #88 of 
| 728.50| 74.49 |-94.06 | 167.08 |-7.56 | |-33.59| 77:69| 538./7|-a/4.53 Not 
| 5:37 |-/4/.09 | | S57 -13 7.05 | 295: 
“CY | 60.96 |-/3.9/ |-27.95 || -—CY | 163.37 |-/78.#/| 130.70 |- 78./0 
$9.82 | (882 |-27.74| 12.0% | -3.47 S |-57 -747 | 23.87) 4/2] 3.0¢ 
one 
S#.9/ |-79.38| 44.43 |-43/ 278.47 | | STATI | MIL 
| /39.27| 93.00 |-/09.24| 26/81 |-3.87 | /755.8/|-2-95 | 95:73 | 
2s -/2.02 | /$¢.29 |-24¢5./3| BST | /38.3f|- 736.92) 376,68 
| 74.70 |-22.72 |-33.4/ | | 2/855 209.0 f|-B0.// ing 
F#-57 | | -15.25| 13.96 | -2a.#0 Ss G46 | 23.77 | 4.9/ | ‘toh 
7 z incli 
| B/.01 |-C-83 | M | 360.23 | | so tk 
83.76 |-75.53 | 39767 |-/46./9 1756.06 | C155 |\525.73|-26 7.98 
-22.0/ | 25.63 |-350.52| 92.60 | | (54.09 32| 340.93 
CY -S5./2 | || -CY [768.08 1266.07 |- 224,76) 30/.80|82,3F The 
& 27.25 | 6.62 |-2.28 | /2.04# | .03 |-/7-80 |-/3.77 | af 
| 220.10 | 96.25 |-V9.7/ | 75:86 | At [272.58 | | 33.// |/78.97 |-y2-22 (For 
| 1000.75 | 54. |-2p.56| 1729.29 |-3-72 | |572.53 278.37 wi 
EIT -33-07 | 233.03) ¢JT ~20f63 | |-98f. 3/ 1373.67 altoge 
-CY |/564.41\ -/8.84 |-=572 | “Cy +/8$/.28 | 306.27 |-23¢.78 | 378.93 |-90.8 No 
five b 


262.29 | /0a-// |-/03.05| /08./2 |-2-B/ | 3779 | M, 


“4 4 
264-96 12.78 | 32.05| 099.45 44.65 |5/2.53 +278.37 


2/7 -#9.93| £6.90 |-556.98) 267.77 | SST -207.3¢| /78.96 |-/020.73| 380.55 But t 
~CY |-1693.25| 29-55 |-64.32| 64.84|-703 | }/8375¢| 3/4.48 |-234.84| $30.50|-94.32 and 1 
6.#7 |-2-4/ | .63| /3.57| 3.3¢ s -2.60 |-1.23 |-77.70| 7-87 must 
eutboard solution tions 
{ 
244-86) $6.65 | 572.53 278.37 
7.36| 178.96 1020-73) 380.55 ut t 
S3C.0/ | -6,50| |-70.53|) 77S It i: 
little { 
M,(0) 
yi(0) = 27.00980691 — 1,947.955958(4, — Y,) = 0 The first and the last four equations readily determine 
(48c) _ the tip slopes and tip deflections gives 
M,(0) = —605,937.1736 (ash) = 0.003169321408, RY; = —0.008342698523R (49) while 
= 0.006012427554, RY, = 0.005946581161R (0) 


y4(0) = —37.39367096 — 3,159.091619%, + 
9,482.335974Y, = 0 (48i) and the numerical columns cos 0, cos y, sin ¥ of Table 4b 
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follow at once, giving the component bending moments 
in lb. ft., the shear components in Ibs. (y and y are of 
minor interest and are not entered). But the first 
harmonic equations result in an inconsistency. 

The explanation of the fact that the Eqs. (48b) and 
(48c) involve only the difference $; — Y, rather than 4, 
and Y; separately, is evident enough. (The arguments 
that follow apply equally to #, and Y2.) Itis seen from 
Fig. 6 that a change of coordinates from yg; and 
(which express slope and deflection with respect to the 
normal plane of the blade) to % and 9; (which express 
slope and deflection with respect to an arbitrarily in- 
clined plane) does not affect the expression for the sum 
of the first harmonic centrifugal and inertia forces. 
Noting that . 


one obtains 


(J: — 6C;z,)/m'dra? = — = — 
(50b) 


and, therefore, with regard to the first harmonic bend- 
ing moment only the difference g¢ — y/r is of signifi- 
cance. In particular, one can choose the arbitrarily 
inclined plane parallel to the tip path plane of the blade, 
so that 


ji, tip = 0 (51a) 
The corresponding tip slope is 
Pr, tip (51b) 


(For the sake of simplicity the bars over y:, Yi, g1, and 
, will be omitted hereafter.) The column Y, is thereby 
altogether eliminated from the analysis. 

Now one can give a clear formulation of the dilemma. 
A solution of the problem is sought which is subject to 
five boundary conditions 


M,(0) = 91(0) = MiQ1) = = = (52) 


But the condition y:(1) = 0 eliminates the Y; column, 
and thus there is only one constant left, #,. This 
must be determined in accordance with the two condi- 
tions 

M,(0) = 0, 91(0) = 0 (6) 


But unless the two expressions 14,(0) and 4,(0) are 
proportional, the problem has no solution. 
It is seen, however, that the two expressions deviate 
little from proportionality. 
M,(0) = 335,653.5109 — 24,223,138.934, = 0 , 
(Yi 0) (53a) 


gives 
&, = 0.01385673062 (53b) 
while 
»:(0) = 27.00980691 — 1,947.9559584, = 0, 
(Y¥; =0) (54a) 


Fic. 6. Reference planes, angles, and deflections for first 
harmonic calculations, P, = normal plane; P, = arbitrarily 
inclined plane; P; = plane parallel to tip-path plane. 


gives 


, = 0.01386571745 (54b) 


Since the moment condition M, (0) = 0 is the cardinal 
root condition, one may use 9; for the calculation of the 
numerical values of ,(x) and S,(x) from the columns 
1, and ®,. The outboard solution so obtained is ex- 
pected to be good at all x except near the origin, where 
it diverges from the true solution to give (0) ¥ 0. 
However, one can produce exact compliance with the 
condition y,(0) = 0, provided a suitable load is applied 
somewhere along the blade. The tip section is probably 
the most convenient place for this purpose, and one 
can call the force thrown into station x = 0.975 as 
Si*. Now there are again three columns on hand: 
1,, ®;, and S,*. (See column at the right end of Table 
4a.) , and S;* are to be so determined that (0) = 
yi(0) = 0 be exactly complied with. The root conditions 


M,(0) = 335,653.5109 — 24,223,138.936,* + 
3,617.3437145;* = 0 (55a) 


yi(0) = 27.00980691 — 1,947.9559586,* + 
0.2908408253.5,;* = 0 (55b) 


yield 
@,* = 0.06081557288, Si* = 314.4546523 (56a) 
and, similarly, 142(0) = y2(0) = 0 yield 
,* = —0.05452102364, S.* = —253.0417395 (56b) 


The inboard solution so obtained is expected to be 
good near the origin, but is bound to diverge near the 
tip from the true solution, since the condition 5,(1) = 
0 (.S2(1) = 0) is violated there. 

The first harmonic results tabulated in Table 4b 
and plotted in Figs. 7 and 8 represent the outboard 
solution from x = 1 to x = 0.6, the inboard solution 
from x = 0.4tox = 0. The tables show both inboard 
and outboard values for x = 0.5;* the graphs show 
their averages. 


* It is to be remembered that the first harmonic 2éJ ani Ce 
are referred to different reference planes, and only their differ- 
ences are comparable. 
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COMPARISON WITH THE COLLOCATION METHOD AND 
WITH JOHNSON-MAyYNE’S QuasI-STaTIC METHOD 


The constant part M(x) and the variable part 
Mi(x) cosy +... M,(x) sin 2p (57) 


of the bending moment M(x) are plotted in Fig. 7a for 
y = 0°, 45°, 90°, ...., 315°, and a special plot vs. y 
is shown in Fig. 7b atx = 0.2. The short dashed curves 
+[(Mi? + (My? + Fig. 7a, consti- 
tute a bound for the variable part which is reached 
only when the maxima of the first and the second har- 
monic bending moments are in phase. In Fig. 8a sepa- 
rate plots are given for Mo, for (Mi? + M,?)'” and for 
(M;? + M,?)'”. Figs. 7a, 8a, and the full line curves 
of Fig. 7b are based on the tabular method. 

The writer repeated the » = O and the »p = 0.35 
bending moment calculations by the collocation method, 
replacing the original m’(x), I(x), c(x) functions by 
smooth distributions, with no extrema or inflection 
points, and having the same Oth, first and second mo- 
ments with respect to the origin as the original distribu- 
tion. (E.g. cotocation = JS, x°dxm' originals etc.) 
For c(x) the value 1.25 ft. was used from c = 0 tox = 
1. (Tip loss was not taken into account.) The de- 
flection curve, y(x), was collocated at the eight stations 
x = 0, 0.1, 0.2, 0.4, 0.6, 0.8, 0.9, 1.0. Curve C of Fig. 
2 represents the collocation bending moment at p = 
0. Because of the absence of tip loss, there is no dip 
below the axis near the tip. A tabular calculation, 
curve D of Fig. 2, performed for the collocation m’, J, 
¢ distributions, shows excellent agreement. (This con- 
tradicts the conclusions reached by the Goodyear group. 
Their poor results for the collocation method must be 
ascribed to a failure to use a station nearer to the tip 
than x = 0.75.) 

The results of the u = 0.35 calculations by the collo- 
cation method are shown in Fig. 8b. Comparison of 
Figs. 8a and 8b indicates that the Oth and the second 
harmonic bending moments are considerably affected 
by tip loss, whereas the first harmonic bending moments 
are only slightly affected. Except for the effect of tip 
loss, the collocation and the tabular calculations are 
again in good agreement.* The phase angles of the 
extrema of the bending moments also agree, as can be 
seen from Fig. 7b which compares bending moment 
vs. y at x = 0.2 by the tabular method and by the 
collocation method. Maximum bending moment oc- 
curs at about y = 0°, a smaller maximum at y = 180°, 
minimum at y = 270°. It is interesting to note in this 
connection, cf. Fig. 9, that maximum flapping of the 
blade occurs near Y = 215°, minimum near y = 30°. 
Fig. 9 is based on the collocation method and shows the 


* In Fig. 8b bending moment curves are shown also for » = 0 
and for w» = 0.175 = 0.35/2, to bring out the fact that M, is 
approximately independent of and M; are approxiniately 
proportional to 4, M3; and M, are approximately proportional to 
u? 


1947 


flapping angle to true scale, deflection of the blade from 
the rigid blade to a 4.6 times enlarged scale. 

The gravity of the second harmonic bending moment 
is rather unexpected. It may be caused by two factors: 
(a) excessive second harmonic air loads, (b) closeness 
to resonance conditions. The author performed a col- 
location calculation, using only first harmonic flapping 
coefficients, but determining bending moments to sec- 
ond harmonics. He found that My and (M,?2 + M,*)* 
remained about unchanged, but (M;? + M,2)'”* in- 
creased by about 50 per cent in the absence of a2 and 
be. So, second harmonic flapping actually effects a re- 
duction in the bending moment. 

The writer was not able to devise a procedure for 
investigating the effect of factor (b). This however is 
readily accomplished by Johnson-Mayne’s quasi-static 
method. Johnson-Mayne, taking cognizance of the fact 
that the maximum bending moment occurs near y = 
0, replace the dynamic air force, Eq. (41), by the ficti- 
tious static air force 


= 6Lo + + TydLs = — 
21 1152] +- =? ul’; (do 1 
— — + + (58) 


where 


1 + 1 


are first and second harmonic amplification factors, and 
v denotes the natural bending frequency of the rotating 
blade as obtained from the relation 


Ex = Ec + Es (60) 


The energy expressions,®: * adapted to the present 
notation, 


tip 
= rx — Hx)*}, 


(61a) 
Eo = = 
Cd(cos rx — H/r)?}, (61b) 
Es = REI (d2y/dr*)2dr = 1/(r/R)*X 


{Eld sin? xx}, 
are based on the assumed deflection curve 
y(x) = sin rx — Hx (62) 


of the blade in the first natural mode, where H has the 
value 


H = sin (63) 


For the McDonnell blade one finds at u = 0.35, 
n = m = 1, that 


t This is indeed the case as can be seen from Table 4b. But 
condition (b) aggravates the situation and is usually the critical 
factor. 
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1 c) moment at = 0.35. Tabular method. 
62) 2.35, N 
ch : Na 
B- QuasisTaTic \ 
53) | C= TiPLossyS, \ 
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But 40 
ical Fic. 8b. Amplitudes of the harmonic components of the bending Fic. 10. Comparison of tabular, collocation, and quasi-static 


moment at » = 0, 0.175, 0.35. Collocation method. bending moments at u = 0.35 for y = 0. 
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H = 1.02878, 2E, = 0.49999»? \ (63) 
2E¢ = 937282/R?, 2Es = 6532800'/R4 


and, therefore, 


v = 63 rad. per sec. = 600 cycles per meal (64) 
lr, = 1.16, Py = 2.22 


and 


= {—0.01381374 — 0.04302451x + 
0.26675537x?}, (65) 


A static tabular calculation, based on 6L,,,, now per- 
mits approximate determination of the maximum bend- 
ing moment of the blade. The Johnson-Mayne curve 
B of Fig. 10 is in excellent agreement with the curve A 
which represents the tabular result of Table 4b for y = 
0. The collocation bending moment curve C is also 
shown for comparison. It is seen that a suitable lower- 
ing of this curve, to account for tip loss, brings also the 
collocation curve into good agreement with the tabular 
curve. 


BENDING MOMENT ANALYSIS OF TEETERED ROTOR 
BLADES 


When the blade is teetered (i.e., when two blades are 
rigidly connected, with see-saw mounting at the hub), 
the calculations are carried out from tip to root by the 
process described in the preceding sections. However, 
at the hub station the slope of the blade cannot change 
abruptly across the mounting. This yields the bound- 
ary condition . 


[¢(0) Jat ver = —[e(O) Jat (66a) 


(which takes place of Eq. (12a) used for individually 
hinged blades) and is to be used in conjunction with 


= 0 (66b) 

Eqs. (66) read in detail 
go(0) = ¢3(0) = ¢s(0) = 0, Mi(0) = M2(0) = O (67a) 
yo(0) = 91(0) = y2(0) = = = 0 (67b) 


(even harmonic slopes and odd harmonic curvatures 
vanish at the origin), and their solution for ®y, ...., Ys 
is no more difficult than that of Eqs. (12a) and (12b). 

It is to be remembered that in using the rigid blade 
air-force approximation for teetered blades, Eq. (44) 
must be replaced by 


M, + Ir(8 + = —MagRec. + 
(44*) 


where the bending moment at the origin 
M, vad Mw + M;s cos 2y + My sin 2y (44**) 


involves only even harmonics. Also, the flapping angle, 
Eq. (42), must contain only odd harmonics, 


=acosy+ (42*) 


As a result, do, d2, b: must be omitted in Eqs. (45) (b 
vanishes as a consequence), and the left-hand sides of 
the equations must be replaced by M,o/R, 0, 0, M,3/R, 
My/R, respectively. 

One might be inclined to seek an extension of John- 
son-Mayne’s quasi-static method to teetered blades by 
using the deflection curve y = cos '/2rx, and making 
corresponding adjustments in the energy expressions, 
Eqs. (61). However, because of the unsimilarity in 
the boundary conditions for even and odd harmonic 
components, the applicability of the quasi-static ap- 
proach to this problem is somewhat questionable. 


ANALYSIS OF THE CHORDWISE BENDING MOMENT OF 
HINGED BLADES WITH DAMPERS 


In determining the chordwise bending moment of a 
blade, Fig. 1 and the earlier calculation schedules again 
apply (M, ¢, ¢’, y, S now refer to bending moment, 
slope, curvature, deflection, and shear force in the plane 
of rotation), but the shear force, cf. Fig. 11, now hasa 
different set of components, to wit: 


6S = 6J — 6C, + 6T + 6C + 6D (68) 
where, in the rigid blade approximation, 
6J = inertia force = —7y/Rw?* (69a) 


—é6C, = lag component of centrifugal force = 
—r{xe — 
6T = lag component of tension force = 
{rx} (69c) 
8@ = Coriolis force = —27xBB/w (69d) 


6D = lag component of air force = 
6D — x6L = '/2pcd(cp — xe,)U? = 
{ (xq + Ur? + — + 
(xo — a)Up?}, (69e) 
Here 
Up/wR = — Bucosy — xB/w (70a) 


Ur/oR'= x + — cosy — (x — 
(70b) 


Fic. 11. Deflection of ee in plane of rotation. Diagram of 
orces. 
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If the flapping angle 8 and the lag angle ¢ of tlie rigid 
blade are known, then determination of the bending 
moment in lag does not differ from that in flapping. 
The calculations are carried out, from tip to the lag 
hinge, by the long or the short schedule. (In the latter 


= — (71) 


replaces 3” { —6C, + 6T}.) At the lag hinge the ap- 
propriate boundary conditions must be imposed. For 
a blade fitted with a viscous damper of D Ib. ft. per 
rad. per sec. damping moment and a centering spring 
of K Ib. ft. per rad. restoring moment, these are (J = 
lag hinge) 


M()) = —D¢(l) — Ke() (72a) 
= 0 (72b) 
The conditions read in detail 


Mi()) = 2wes(l) — Kex(2) (78a) 


= .... = = 0 (73b) 


Their solution for yp, ...., Y4 proceeds in the same 
way as for the Eqs. (12a) and (12b). 
Determination of the lag angle 


Aicosy + sin y + Az cos 2y + sin 
(74) 


is effected from the equation of moment equilibrium 
with respect to the lag hinge, 


+ we] = —201,(1 + )68+Q, (75) 


The flapping angle 8 which enters both the expression 
of the Coriolis force and the expression of the lag torque 
Q, is a quantity that is already known from the solu- 
tion of the flapping problem. (Strictly speaking, the 
¢and ¢ terms in Uy, Eq. (70b), modify somewhat the 
flapping coefficients, but this effect can be neglected.) 
Eq. (75) decomposes into five rather cumbersome com- 
ponent equations from which Ao, Ai, ...., Bz can be 
determined. 

The author has not as yet carried out bending mo- 
ment calculations by the tabular method. However, 
by the more cumbersome 8-station collocation method 
results are available for the 23-ft. blade at uw = 0.35 
for D = 400 Ib. ft. per rad. per sec., K = 0, and for D = 
K = 0. Since disregard of tip loss has no counterpart 
in the plane of rotation, the collocation method is ex- 
pected to give reliable results for the chordwise bending 
moment. 

Neglecting the ¢ and ¢ terms in Eq. (70b) and in 
Q, (the error so committed is thought to be moderate), 
and using 


R, = 1.3958 ft., I, = 328.4 slug-ft.’, el, = 
34.36 slug-ft.2 (76) 


Eq. (75) leads, for D = 400 Ib. ft. per rad. per sec., 
to 


Ayo = 0.149964, Ai = —0.021006, B, = 0.028314, 
Az = 0.004669, B, = —0.001806 (77) 


and for D = 0, to 


Ao = 0.149964, A; = —0.022656, B, = 0.027090, 
A, = 0.004718, Bz = —0.001681 (78) 


The differential equation of blade deflection, ob- 
tained from 


6°M/ér? = 6S/ér (79) 


assumes the form 


dr? dr? 


(? - + m'y = —2m'wBBr + (80) 
dr r dr 


2 
(e122) — 1?) + mira? X 
r 


After replacing the variables r, y in Eq. (80) by 
u = (r— R,)/(R — R,) (81a) 
T = y/R + Ticosy+.... + sin 2y (81b) 
and defining the polynominals 


P, = — + a, = RY Di. 
(82b) 


the deflection functions 


T = anP.(u) + Cy"Pi(u) + Cs"P3(u) + 
CrP(u) + .... + Co"Po(u) (83) 


are collocated at the eight stations, u = 0, 0.1, 0.2, 0.4, 
0.6, 0.8, 0.9, 1.0. This fixes the values of the 40 coef- 
ficients C,°, ...., 


The result of the calculations, the chordwise bending 
moment M(x) is Shown in Fig. 12a for y = 0°, 45°, 


T 
bd 
| 
40 


Fic. 12a. Bending moment in plane of rotation with hinge 
damper, at » = 0.35, for various azimuths. 
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R, = So/ Fo; R, = fi/ Fi (85) 
Ty | | the two stress ratios, the margin of safety 
MS. = (Ro + 1 (86) 
| | ral must be a positive quantity when My, and M, are the 


Fic. 12b, Bending moment in plane of rotation without damper, 
at » = 0.35, for various azimuths. 


...+, 315° when D = 400, and in Fig. 12b when D = 
0. It is seen, cf. Fig. 13, that the blade with damper 
has a maximum hinge bending moment of about 400 Ib. 
ft., and that the bending moment rapidly reduces, as 
one moves outboard along the blade, to the moment of 
the blade without hinge damper. ‘ 

It is to be noted that when D is small compared to the 
critical damping for lag motion, the deflections ¢ for 
D > 0 and for D = 0, Egs. (77) and (78), are nearly 
equal. It follows that for small values of D the lag 
hinge bending moment M(/) is approximately propor- 
tional to D. 

STRESS ANALYSIS 


The blade spar is subject to spanwise bending 
stresses, chordwise bending stresses, and centrifugal 


tension. 
The critically loaded fiber (at station x) is thus sub- 
ject to a combined stress of 


folx) + fil) = [Mr%(x) + + 
c(x)/A(x) (84) 


With F, denoting the ultimate allowable bending stress, 
F, the ultimate allowable tensile stress,’ 


at 
Dad \ 
\ \ / 
. 
atuso 360° 


-200 
with damper 


ith out damper 
1 l 


Fic. 13. Bending moment in piane of rotation at stations u = 0 
and u = 0.2 with and without damper. 


at rated w. 


maximum bending moments in flapping and lag calcu- 
lated for rated w, 2g load, and multiplied by a factor 
of 1.5 to give 3g limit bending moment, and by an 
additional factor of safety of 1.5 to give ultimate 
bending moment;* and when C is calculated for 50 
per cent overspeed and is multiplied by a factor of 
safety of 1.5. 

In checking the spar for fatigue, 1g stresses are used 
Although the two extreme stresses 


maz. fi + So, mar.» fi + So, min. (87a) 


(f is positive for tension, negative for compression) are 
not reached for the same fiber, it is conservative to dis- 
regard this fact, and determine the margin of safety in 
fatigue by the formula 


where 
Su = + fmin.)» Sa = f maz. —Sn (87b) 


are the mean and the alternating stresses, respectively; 
K is the stress concentration factor, and F,, and F, are 
the allowable steady and alternating stresses, deter- 
mined in accordance with reference 10. 

The above procedures outline the principal features 
of the strength check method which was used by the 
author’s colleague, A. A. Lischer, in the stress analysis 
of the McDonnell 23-ft. blade for flight attitudes. 
Satisfactory positive margins and were 
found for all cross sections of the spar. 

It is now felt that the specification in static analysis 
of a 3g limit load factor for f,, 50 per cent overspeed for 
fs is overly conservative. Recent investigations indi- 
cate that the assumption of a maximum overspeed of 
25 per cent is more in line with facts and that the maxi- 
mum maneuvering load factor of the McDonnell heli- 
copter cannot exceed m = 2.6 and the maximum gust 
load factor cannot exceed 1.9. The latter results were 
obtained by the author’s colleague, A. Gail, in an exhaus- 
tive (unpublished) study of the load factors to which 
helicopters may be subjected. 
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* A direct calculation for 3g load-is not permissible because the 
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3 for limit bending moment and the additional safety factor of 1.5 
give unduly conservative results, cf. Fig. 5a. 
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Letters to the Editor 


Dear Sir: 

I have read with great interest the three papers by Georges 
Welter on the buckling of curved plates published in the July, 
1945, and April and November, 1946, issues of the JouRNAL. 
I should like to congratulate the author on a most important 
technical advance and reproach him for an almost equally impor- 
tant scientific regression. 

Following the work of von Karman and Tsien,! of Leggett,? 
and of Koiter,* the elements of the theory of the buckling of 
curved plates are no longer open to debate; the general outline 
may be said to be internationally agreed, and the few details 
that remain in dispute do not affect the question that the author 
attempts once more to put at issue. It is indeed the most satis- 
fying aspect of the theory that it clearly indicates the vital im- 
portance of initial irregularities and thus accounts for the ex- 
traordinary variability of behaviour of curved plates which has 
worried experimenters for years. In order to coordinate and 
interpret the author’s own results, reference to the initial ir- 
regularities cannot be avoided, and indeed I would claim that the 
chief value of the author’s work lies in the means he has devised 
to control this factor. In this respect the author has succeeded 
where the writer and several others have failed; but the very 
basis of his success rests in just those initial irregularities that he 
dismisses so cavalierly. 

The detail of the analysis of the postbuckling behaviour of a 
curved plate is extremely complicated so that the full treatments 
of references 1, 2, and 3 arefarfromeasytoread. Asaresult the 
elements of the theory are perhaps not yet widely understood, 
and to discuss the effect of initial irregularities it is necessary 
first to rehearse the theory in outline. 

Suppose that in a test on a flat plate under compression we 
plot depth of buckle W (measured at any arbitrary point) against 
compressive strain e; we shall find that W? ae — e,. Some excel- 
lent demonstrations of this relation have been obtained during 
tests made by D. J. Farrar and his collaborators at the Bristol 
Aeroplane Company. . 

By reference to theory also it is quite easy to show that so 
long as the wave form of buckle remains unchanged, the curve 
BAC in Fig. 1 is a true parabola symmetrical about the axis of 
strain. Ignoring the singularity that a perfectly flat plate should 
never buckle at all, we may say that in a perfect plate inward 
and outward buckles are equally probable; whereas any practical 
plate must have slight initial irregularities that predispose it to 
buckle one way or the other. As a result for any practical plate 
the W — e diagram follows some line such as ODC instead of 
OAC. 

To each and every conceivable buckled wave form there corre- 
sponds a parabola like BAC in Fig. 1, and the true buckled form is 
that which corresponds to the parabola that intersects the axis 
of strain at the lowest value (Fig. 2). Nevertheless, if the plate 
has initially only a small amplitude of buckle in the mode corre- 
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sponding to the curve BAC and a much larger amplitude of buckle 
in a mode corresponding to a curve B’A’C’, where OA’ is slightly 
greater than OA, the mode B’A’C’ may be developed in prefer- 
ence to the mode BAC. This point again has been well demon- 
strated in the tests made at the Bristol Aeroplane Company al- 
ready cited. Transition from one wave form to another after the 
final form is well developed is occasionally possible,‘ and this proc- 
ess is similar to the buckling of initially curved plates. But it 
should not be confused with modification of the wave forms by 
gradual development of harmonic terms. 

When the plate is initially curved, the chief and almost the 
sole effect’ on the W — e relation is to shift the parabola to the 
left (Fig. 3). Then, according to the sign of the initial amplitude 
of buckle in the mode corresponding to the curve BAC, either 
ODB or OEC is possible. The great practical importance of this 
difference is indicated in Fig. 4, which shows the stress-strain 
curves corresponding to the two branches XC and XAB of the 
wave depth curve. (Incidentally the curve CXAB in Fig. 4 is 


also a true parabola so long as the wave form of buckle remains 
unchanged and it touches the elastic line OX at X. The width of 
the parabola in Fig. 4 is proportional to the amount by which the 
parabola in Fig. 3 is offset from the centre line, and both are 
determined by the value of b?/Rt, where b is the width of the 


1947 


plate, ¢ its thickness, and R the radius of transverse curvature, 
For a flat plate the parabola in Fig. 4 is indefinitely narrow and 
the two branches XC and XAB become indistinguishable.) 

Unfortunately, “of all the possible buckled forms there is a 
wide variety, which, when represented by curves such as BAXC 
in Fig. 3, intersect the axis of strain close to the minimum value.§-8 
Accordingly, in order to suppress the branch XAB in Fig. 4, it is 
necessary to secure a positive value of the initial amplitude not 
merely of one wave form but of several. 

Previously, this desired end has been sought by making the 
plate slightly barrel shaped® but without success. It is clear from 
the stress-strain diagrams obtained by the author that he has 
achieved the desired end, and it is understandable that the process 
of overforming could have this result. Unfortunately, this 
explanation scarcely suffices to cover all the results obtained by 
the author. 

The chief difficulty is that the critical stresses for primary 
buckling recorded by the author are nearly all much higher than 
the values obtained by other experimenters. The results re- 
ported in R. & M. 1894 (to which the author refers in his first 
paper) were made on plates that were not well formed, and the 
values of critical stress obtained may be expected to be rather 
lower than the average. On the other hand, the values obtained 
in recent tests made by K. B. Jackson and A. H. Hall at Ottawa 
for the National Research Council of Canada” are still well below 
the theoretical values, and these plates and the tests on them were 
made with exceptional care. 

From the general discussion above it is clear that any inward 
buckle causing a fall of load at buckling should occur well below 
the theoretical limiting value; this value is 0.6E(t/R) for a com- 
plete tube or wide panel. Although in a narrow panel the effect 
of restraint at the sides raises the critical stress, at the values of 
b?/Rt of about 40 and over used by the author this effect is prob- 
ably not appreciable. This is indicated by Redshaw’s curve 
shown in Fig. 19 of the author’s first paper. 

It is conceivable that the edge supporting members used by 
the author may have carried an appreciable proportion of the 
total end load, but this possibility seems to be ruled out by the 
tests made with wide clearance in the slots (Fig. 8 of the second 
paper). The only other explanation that suggests itself is that the 
wave form of buckle was restricted in some special way. The 
length of the plates tested was less than that used in other tests, 
but this difference again seems to be covered by the results re- 
ported in Fig. 11 of the first paper. 

In view of the technical importance of the overcurving process 
described in the third paper, it seems most desirable to explain 
the discrepancies between the author’s results and other work. 
Could the author carry out some tests on complete tubes or much 
wider plates, so as completely to eliminate any possible restric- 
tion on the wave form due to the side supports? 


H. L. Cox 
Engineering Division 
National Physical Laboratory 
Teddington, Middlesex, England 
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Dear Sir: 

It was with deepest interest that I read the discussion of my 
three papers concerning ‘‘The Buckling of Curved Aluminum Al- 
loy Sheets,” published in the JouRNAL in 1945-1946. I am 
thankful to Mr. Cox for taking up the question of the experimen- 
tal results, recorded at Ecole Polytechnique, with regard to the 
existing theory on the buckling phenomenon of curved plates. 
Iam also much obliged to Mr. Cox for his authoritative statement 
that a most important technical advance has been achieved in 
this field. Furthermore, Mr. Cox must be commended for his 
outline of the elements of the theory presented for a better under- 
standing of the effect of initial irregularities on the buckling loads 
of thin plates. 

In view of the fact that the experimental results of this inves- 
tigation have met with some interest in the circles concerned, 
a few words about the historical development of this work, car- 
ried out at the beginning of the war with limited funds and in a 
rather short period of time, might be of help to better understand 
the scope, as well as the limitations, of this investigation, the 
published report of which was condensed by more than 50 per 
cent of its original form. 

The experimental part of this investigation was begun at the 
end of April, 1942, and was completed by the beginning of No- 
vember, 1942. From this time-period we must deduct about 6 
to 8 weeks, during which the universal Baldwin testing machine 
was given over to another equally important work. This means 
that the main experimental tests of this investigation were carried 
out in less than half a year. About 20 to 30 per cent of this time 
was occupied with preliminary work before starting the main test 


series. These preliminaries consisted in a complete verification 


of the accuracy of the testing machine; a check of the transla- 
tional movement of the tables under load; the determination of 
the stiffness of the machine; the adaptation of a special high- 
precision load-strain recorder; the development of the whole 
testing technique, such as the gripping of the sheet specimens, as 
well as the development of different kinds of restrainers; the 
construction of new devices; the preparation of special speci- 
mens, etc. Because of the highly urgent character of this prob- 
lem, the experimental work had to be accomplished during the re- 
maining short period of time. Also, within the same time limit, 
other basic work on an equally important subject had to be in- 
cluded. This concerned the investigation of the directional 
characteristics of the material tested under tension and com- 
pression loads, such as Young’s modulus, elastic limit, micro- 
deformations, the effect of Bauschinger with single thin test speci- 
mens, internal stress analysis of residual stresses in rods and 
sheets, etc. The results of this investigation were published in 
the Proceedings of A.S.T.M., 1944, and this paper was one of 
the first reports in this important new field. These tests were 
necessary for a better understanding of the behavior of available 
sheet material in panels under buckling loads. It is evident that 


under these conditions a simple yet effective test procedure had 
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to be developed without sacrificing any of the required precision. 
The experimental results are founded exclusively on the auto- 
matic, high-precision records of the testing machine itself, the 
test procedure having been systematically checked for any pos- 
sible errors. The main experiments were carried out personally 
by the writer, with the assistance of one or two men only. The 
aim was to develop an independent automatic recording of results 
in the shortest possible time under restrictive economic condi- 
tions. This was achieved to a satisfactory degree. By this 
means, the recorded results are, as far as possible, free from the 
influence of any variable human factor, and the only way to re- 
port upon these results in a most objective way was to show the 
original diagrams recorded by the Baldwin machine, as was done 
in the three publications here referred to. In face of these facts, 
it is hard to see where the point of a scientific regression, as 
charged by Mr. Cox, may come in. This regression, if any, is 
completely beyond my control, and this criticism may eventually 
be charged against the Baldwin machine itself, the machine 
which recorded such unprecedented high buckling loads. If the 
results recorded in these three papers can be ignored, then the 
question of any regression can be equally disregarded. If, how- 
ever, these results will prove to be consistent and if, by further 
tests, no secondary effects such as the influence of supporting mem- 
bers and the dimensions of the panels, etc., can be found by which 
indirect means the curved plate may be strengthened, then we 
have to consider (and this will not be the first time) how the for- 
mula can be changed to fit the results. This will have to be done 
even though, as stated in my first report, still better results than 
those already shown are attained by further development and 
improvement of the test equipment. What now seems to be a 
scientific regression may, in the light of higher results, be revealed 
as an important scientific advance. Unfortunately, despite the 
encouraging results obtained, these tests had to be abandoned at 
the end of 1942 because of lack of funds, and any investigation re- 
garding the influence of initial irregularities had to be, as Mr. Cox 
puts it, cavalierly dismissed. 

For further investigation of this kind special funds are neces- 
sary. This further research would include some fundamental 
work applied to the initial irregularities of the panels, which 
might possibly contribute to the observed higher buckling loads. 
To find an adequate answer to this special problem may absorb 
more time than the main work itself. I am in complete accord 
with the statement by Mr. Cox that, in view of the technical 
importance of the overcurving process, it seems most desirable to 
explain the discrepancies between my results and those of others. 
I agree also to resume some of these experiments, under such con- 
ditions, however, which would seem to be favorable to a higher 
buckling strength in order to find the main cause for these unex- 
pected favorable results, which, in fact, may be of vital technical 
importance relative to the behavior of curved plates in stressed- 
skin constructions of modern aircraft. 

GEORGES WELTER 
Professor of Applied Mechanics 
Ecole Polytechnique de Montréal 


Dear Sir: 

In his paper on ‘“‘Aerodynamic Characteristics of Rectangular 
Wings at Supersonic Speeds” in the February, 1947, issue of the 
JourRNAL, E. A. Bonney did not seem to have considered the re- 
duction of wave drag at zero angle of attack—i.e., the ‘form 
drag’’—at very small aspect ratios. Therefore, when his results 
are applied to design problems, one must be aware of this limi- 
tation. 

That there is such a reduction in form drag can easily be seen 
for the case of untapered wing by following the analysis of M.. J 
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Lighthill (R. & M. No. 2001, ARC). Let c be the chord; 2b, 
the span; and a = \/M?— 1. Then there is no overall reduc- 
tion in form drag if c/a < 2b, as shown by Lighthill. However, 
if c/a > 2b—i.e., when the Mach cone from the right tip of the 
leading edge includes the left tip of the trailing edge and vice 
versa—then the drag reduction can be calculated as 


ab b 


— (x/a) 
fe) dé 
— (x — — a? (b — y)? 


This formula is correct if the span is much larger than the thick- 
ness of the airfoil. f(x) is the source distribution function given 


by 


f(x) = —(U/2x) (dt/dx) 
where U is the free stream velocity and ¢ is the thickness of the 
wing at distance x from the leading edge. By an interchange of 
the sequence of integration, 


ab ab 2 
ap = dx dt = at 
b 0 


If Cpr, is the form drag coefficient of a wing of infinite span 
and if Cpr is the form drag coefficient of a wing of aspect ratio 
R = 2b/c, 


In case of a symmetrical parabolic airfoil, 
= 46(x/c)[1 — (x/c)] 

where é is the thickness ratio of the airfoil. For a given 5 and a 
there is a value of R for which Cpr is minimum. This optimum 
R is given by aR = '/; or 2ab = c/3. The corresponding form 
drag coefficient is 

CoPopt. = (16/3) (5?/a) [1 — (4/9x)] 
or a reduction of 14 per cent. 
H. S. Tsren 


Guggenheim Aeronautical Laboratory 
Massachusetts Institute of Technology 


Dear Sir: 

The paper ‘‘Exact and Approximate Solutions in Two-Dimen- 
sional Oblique Shock Flow,” by E. V. Laitone, which appeared 
in the JOURNAL OF THE AERONAUTICAL ScrENcEs for January, 
1947, contained one expression that appears to be incorrect. This 
suspicion has since been confirmed by Mr. Laitone. We therefore 
submit the following comments. 

Laitone pointed out that Busemann’s approximate solution for 
oblique shock flow is incorrect in the third term. It seems in 
addition, that Busemann’s approximate solution for isentropic flow 
is also incorrect in the third term. We believe that the correct 
form of the coefficient that Laitone designated as C; should be: 


1 1+ — —5 ‘ 
CG; (Me? mal 6 Me + 6 + 


5 4 
M, 2] 2 = 
+ 7) M..* + 
subtracting 53; (given by Laitone) from the correct C;, we get 


y+1\ Ma! 5 — 37 6 — 2y\? 
p= ( 12 4 - a 


(5 — 3y)? 


where D is the difference between the third terms for the isen- 
tropic pressure change and the oblique shock pressure change— 
ie., D = Cy — ds. 

These correct values of C; and D at various M should be as 
follows: 


M ee a 1.3 1.4 1.5 1.6 
C; 867 53.8 14.40 5.80 3.06 1.97 
D 322.6 0.58 —0.13 —0.328 —0.271 —0.183 

M 2.0 3.0 4.0 5.0 10.0 

C; 0.927 1.13 1.51 1.91 3.9 

D —0.092 0.061 0.105 0.136 0.322 


This corrected value of C; will bring the approximate expres- 
sions for pressure rise in closer agreement with an exact computa- 
tion. 

ALLEN E. Puckett 

Tine Yi Li 
Daniel Guggenheim Aeronautical Laboratory 
California Institute of Technology 
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Numerical Methods for the Calculation of 
Elastic Instability 


BRUNO A. BOLEY* 
Polytechnic Institute of Brooklyn 


SUMMARY 


Three numerical methods for the evaluation of buckling 
loads requiring an operations table similar to that used in South- 
well’s relaxation procedure are developed and described: the 
Determinant Method, the Convergence Method, and the Energy 
Method. The Determinant Method makes use of the fact that 
the determinant made up by the elements of the operations 
table vanishes at buckling. The Convergence Method states 
that a systematic relaxation process based on the operations 
table is convergent only if the structure is stable. The principle 
that the second variation of the total potential is zero during 
buckling forms the basis of the Energy Method. The theoretical 
background of the three methods is discussed after a review of the 
setting up of operations tables and of Southwell’s method in 
general. The procedure to be followed in each method is given 
and illustrated in numerical examples. Experiments on sheet 
and stringer combinations carried out at the Polytechnic In- 
stitute of Brooklyn Aeronautical Laboratories (PIBAL) are 
described, the results of which are in good agreement with those 
obtained from each of the three methods of calculation. 


INTRODUCTION 


f bes BUCKLING OF A STRAIGHT BAR under axial 
compressive loads is probably the best-known 
examplé of elastic instability. The critical load is 
found in this case directly from the characteristic value 
of the governing differential equation of the problem. 
In more complex structures, however, it is not con- 
venient or possible to set up such a differential equa- 
tion; then strain energy methods in which an analytical 
expression for the buckling distorted shape is assumed 
may be employed to good advantage. In several 
instances, application of such strain energy methods 
may in turn become too difficult or too long. In this 
paper three numerical methods are described which 
may be advantageous in the solution of complex 
buckling problems. 

A knowledge of the basis of Southwell’s relaxation 
procedure!~* is essential to the understanding of these 
numerical methods. Accordingly, the first section 
of this paper contains a review of Southwell’s method. 
The theoretical background of the methods is dis- 


Presented at the Structures Session, Fifteenth Annual Meeting, 
I.A.S., New York, January 28-30, 1947. 

* Instructor, Department of Aeronautical Engineering and 
Applied Mechanics. Paper based on a thesis submitted in 
partial fulfillment of the requirements for the degree of Doctor 
of Aeronautical Engineering at the Polytechnic Institute of 
Brooklyn. The author is indebted to Dr. N. J. Hoff, his thesis 
adviser, for the criticism and advice received. 
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cussed in the next section, while the procedures to be 
followed in each method are given in a third section. 
Results of an experimental investigation are then given 
and are compared with results of numerical examples 
in which the three theoretical methods are applied to 
the experimental specimens. 


SOUTHWELL’S RELAXATION METHOD 


Assume the structure under consideration to be 
made up of several “units” of finite dimensions. 
Typical portions of the structure should be chosen as 
units, such as the bars of a framework or a panel of 
sheet with its edge-reinforcing elements in a sheet and 
stringer combination. Consider now some well-defined 
boundary points, known as “‘joints,’’ in these units 
(for instance, the ends of the framework bars or the 
corners of the sheet panels). If lines are drawn be- 
tween these points wherever direct structural con- 
nections exist, a more or less regular lattice is formed. 
Assume that for simplicity’s sake a particular example 
is chosen which contains only six joints, numbered 
and arranged as shown in Fig. 1. 

All joints but one will, at any one moment, be as- 
sumed to be rigidly fixed in space, and the effect of 
displacements of the one movable point will be in- 
vestigated. This investigation, discussed later in 
more detail as part of the so-called ‘unit problem,” 
must be done for all points in the structure—that is, 
until every joint in turn has been considered movable 
with all others fixed. Let the generalized force exerted 
on joint 2 by a unit generalized displacement of joint 1 
in given directions be denoted by an, that exerted on 
joint 1 by a unit displacement of joint 2 by ay, and 
similarly for all the others. In a more general case, 
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Fic. 1. Example of structure for setting up an operations table. 
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displacements in different directions and rotations 
(and therefore also the corresponding forces and 
moments) should be considered separately. That is, 
ds, may denote ‘the force in the x-direction caused at 
joint 5 by a displacement in the x-direction of joint 6, 
while a;; may denote the force in the x-direction caused 
at joint 5 by a rotation at the same joint 6. 

Referring again to the example of Fig. 1, it should be 
noted that as a consequence of the assumption of rigid 
fixation of all points, with the exception of the one 
which is at the time being moved, displacements of 
points 5 and 6, say, will have no effect on points 1 and 
2 (or a5 = Aig = G25 = do = 0). Assuming the validity 
of the principle of superpositions, the total force 
exerted on joints 1 and 2, for example, after all joints 
have been moved, may be written as: 


aux: + + + =0 (1) 
+ + + + P=O0 (2) 


joint 1: 
joint 2: 
where the x’s represent the displacements of the 
various points of the structure, P is a known applied 
load, and the equality to zero is the requirement for 
equilibrium of joints 1 and 2. Similar equations may 
be written for every point in the structure, and there 
will then be as many equations as there are displace- 
ments. A set of linear equations is thus formed, with 
displacements as unknowns, which represents the 
equilibrium conditions for the whole structure. A 
typical such equation for the ith joint may be written as 


+2 = 0 (3) 


where n is the number of points in the structure. 

The array, or matrix, formed by the coefficients of 
these equations is called the “‘operations table” and 
may be represented in a general case by 


a2 Qo2..... don 
an Ag..... Qin 


It should be noted that, as a consequence of Max- 
well’s reciprocal theorem, = 

The individual elements of the operations table are 
obtained from the solution of what is called the ‘unit 
problem.’ The unit problem is the evaluation, within 


a unit, of the forces and moments in the various direc- 
tions and planes, respectively, caused at all joints by 
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To is the applied torque vector. 


Fic. 2. Curved bar in torsion. 


unit displacements and rotations of each point sepa- 
rately. These forces and moments are known as the 
“<nfluence coefficients,’’ and the elements of the opera- 
tions table themselves are then simply sums of influence 
coefficients. 

The unit problem may be evaluated by any con- 
venient method. ‘As an example, the influence coeffi- 
cient is obtained here, with the aid of Castigliano’s 
theorem, for the case when the unit is a segment of a 
curved bar under torsion. Referring to Fig. 2, the 
moment 7 required at any section for equilibrium is 
T = Tp», where 7) is the applied torque. Considering 
separately the components of T causing torsion and 
bending, the total strain energy U may be written as" 


U = (1/2) fo® [(T cos 6)?/GC]rdo + 
(1/2) [(T sin (5) 


The total angle of twist w = dU/dT>) by Castigliano’s 
theorem; then, assuming the radius 7, the torsional 
rigidity GC and the bending rigidity EJ to be constant, 
the influence coefficient is 


/ + 28) + 


(§ sin 28)] ©) 


Use of Eq. (6) is made in the operations tables appearing 
in the last section of the paper. 

Influence coefficients for beam columns, which will 
be ‘needed in the numerical examples, are assembled 
in Table 1. They were obtained from the differential 
equation of the curve of deflection of a beam in bending, 
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TABLE 1 


BEAM- COLUMN INFLUENCE COEFFICIENTS 


FORCES 
DISPLAC.| UNIT Fy M, 
END A -P -P P 
0, 2 (1-c) sku S = SIN KL 
1 - (kif) = coskt 
k*= t= TaANkL 


SIGN CONVENTION 
FORCES ON CONSTRAINTS 


x L 
dF 


P POSITIVE AS SHOWN 


in which the moment term equgl to the product of the 
axial end load P and the deflection was considered. It 
may be said that in the determination of stress dis- 
tribution by Southwell’s method such terms may 
usually be neglected (see numerical examples), while 
it is necessary to include them for the evaluation of 
stability problems. 

A general discussion of influence coefficients may 
be found in reference 1, which also contains specific 
examples related to straight bars and beams with and 
without axial end loads. Beam-column influence 
coefficients are also given in reference 4. Flat panels 
of thin sheet are treated in reference 3, curved bars 
and curved thin sheet panels in references 5 and 6. 
A detailed discussion of the setting up of operations 
tables and the Southwell method in general is given in 
references 1 and 3. 

Evaluation of the stress distribution in a given 
structure requires the solution of the linear equations 
represented by the operations table. This may be 
obtained by an exact method as, for example, that of 
reference 7, or by an approximate “relaxation’’ pro- 
cedure such as the one shown in the numerical example 
for the Convergence Method. The displacements of 
each point found by such a solution may then be used 
to calculate the stresses in all the members of the struc- 
ture as illustrated in the last section of this paper. 

The ‘theoretical background of methods for the 
evaluation of buckling loads, in which use is made of 


an operations table, is discussed in the next section. 
The actual procedures to be followed in three such 
methods are given after the theoretical considera- 
tions. 


THEORETICAL CONSIDERATIONS 


Let A be the square symmetrical (a;; = a,,) matrix 
shown in Eq. (4), the order of which is . Let r be its 
rank—that is; the number of linearly independent 
rows (or columns). If r = m the matrix is said to be 
nonsingular, if r < m it is said to be singular. 

Let Q be a quadratic form whose coefficients con- 
stitute the matrix A; then 


Q = + + + 
+ + ..... + + 
(7) 
or 
Q=> (8) 


A quadratic form Q is said to be positive definite if 
it is negative for no values of the variables x, to be 
negative definite if it is positive for no values of the 
variables x, and to be indefinite if it varies in sign. 
Some of the most useful properties of quadratic forms 
are collected in Table 2. The term principal minor 
A, found in this table refers to the 7-rowed determinant 
which is a minor of the matrix A, whose main diagonal 
coincides with part of the main diagonal of the given 
matrix A and whose elements are taken from the 
first t rows and columns of the given matrix A. Thus 


(9) 


Additional information about matrices and quadratic 
forms may be found in reference 8. 

The statements contained in Table 2 referring to 
states of equilibrium to which a quadratic form corre- 
sponds are proved in the theorems that follow. The 
elements of the matrix A may be assumed identical 
with those of an operations table corresponding to a 
physical structure. 


Theorem 1 


The matrix, say A, of the elements of an operations 
table is that of a negative definite nonsingular quadratic 
form if, and only if, the structure it represents is stable, 
provided the elements express (as is usual) forces on 
the constraints. 
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TABLE 2 


n n 
Properties of Quadratic Forms Q = > > QijXix;; (aij = aji) 


i=1j=1 
Negative Positive 
Definite Indefinite Definite 
Form Form Form 


Q Nonsingular (Then r = n, A, # 0) 


Q<0 For all x; # 0 For some values of Never 
x 
Q=0 If and only if x; not necessarily If type id only if 
zero 
6 Never some values of For ‘all xi #0 
Principal Positive if 7 is May be positive, Always posi- 
minors A; even, nega- mnegativeorzero; tive 
tive if i is but if A, = 0, 
odd the rows (and 
(Ao = 1 by columns to pre- 
definition) serve symmetry) 


may be _ rear- 
ranged so that 
<0 


Main diagonal All negative May be positive, All positive 


elements negative or zero 
= A;) 
Corresponds to Stable Unstable Unstable 
equilibrium 
Q Singular (Then r < n, A, = 0) 
Q<0 For some Forsome values of Never 
values of x; =X 
Q=0 x; not neces- x; not necessarily x; not neces- 
sarily zero zero ily zero 
Q>0 Never For some 
values of x; 
Principal One A, ¥ Oat least can be found; all A; = 0 
minors A; ifi >r. Fori < the rules hold as given 


for Q nonsingular 
Corresponds to Neutral Unstable 
equilibrium 


Unstable 


An expression for the total potential energy in the 
structure is required for the proof of this theorem and 
is therefore derived here. 

If x; represents the total displacement of point 7 in 
a given direction due to the application of external 
loads P to the structure under consideration, the 
total work Wz, done by these loads may be written as 


n 
We = >» Pix; (10) 
where is the number of points in the structure. By 
definition, the negative of Wg is the potential of the 
external forces. 
The force acting internally between the ith joint, or 
constraint, and an element of the structure at that 


point may be represented in terms of the influence 
coefficients a,, by the quantity }> a,x; The internal 
j=l 


work done at point 7 by this force during the displace- 
ments x; caused by the application of the external 
loads P is 


x n 1 n 
Wi -f ( ay, 7% D + 
fat 


(11) 


= j in the next to the last summation, and 


where 7 


1947 


7 ~ j in the last. The total internal work W, may 
now be obtained by summing up over all the joints. 


Then 


t=1 j=l 
Gi (i¥7) 


where the factor of !/, has been introduced in front 
of the last term because the double summation in that 
term includes twice each element of the structure. 
In Eq. (12) again terms in which 1 = j have been 
kept separate from those in which 1 ¥ j; however, 
as the multiplying factor is now the same for all terms, 
the total internal work may be written as 


By definition, the negative of the internal work W, 
is the total strain energy in the structure. 
The total potential V is the sum of the potential of 
the external forces and the total strain energy and may 
therefore be written as 


= —(Wz+ W,) = Ms > (14) 
i=l j=l 


The well-known relations between the second variation 
of the total potential and the state of equilibrium may 
now be used in the proof of theorem 1. 

The change in total potential due to a small variation 
in the displacements is 


6V = > > aus) bx, (15) 
i=1 j=l 


Incidentally, the quantity 6V is zero when the structure 
is in equilibrium, as can be noticed upon comparison 
with the equilibrium conditions exemplified by Eq. 
(3). If the displacements are now again varied while 
the external loads P; are kept constant and some 
nonzero variations of the displacements can be found 
which again correspond to zero variation of potential, 
but none can be found which decrease the variation 
of potential, the equilibrium is neutral. If all the dis- 
placemént patterns not identically equal to zero in- 
crease the variation of potential, the system is stable; 
if some can be found which make this variation smaller, 
the system is unstable. Hence, for stability for any 
arbitrary nonzero variation of displacements, 


(16) 


Making use of the equilibrium conditions represented 
by Eq. (3), a necessary and sufficient condition for 
stability is found to be 


de, <0 an 
i=l j=l 


— 
for 
forn 
Simi 
Ci 
tabl 
defit 
by 
conc 
if th 
the 
is pc 
| is th 
tabl 
Ce 
if, t 
elem 
crite 
cons 
nece 
one 
repri 
than 
than 
mus 
Si 
nece 
a stipt 
next 
| of tk 
of p 
fact 
| Thec 
Le 
latio 
| for a 
‘ n n n the | 
i=1 j=l j=l then 
| n Ther 
| scrib 


May 
joints, 


7) 


CALCULATION OF ELASTIC INSTABILITY 341 


for all values of the variations not identically equal to 
zero. Therefore, the quadratic form represented by 
the left-hand side of inequality (17) is nonsingular 
negative definite. As its matrix is identical with that 
formed by the operations table, theorem 1 is proved. 
Similar proofs are given by von Mises® and Hoff.’ '! 

Corollary 1.—-For unstable equilibrium the operations 
table corresponds to either an indefinite or a positive 
definite quadratic form. 

Corollary 2.—The value of the determinant formed 


by the elements of an operations table under stable 


conditions is negative if the order is odd and positive 
if the order is even. For an unstable system between 
the first and second buckling loads, the determinant 
is positive if odd-ordered and negative if even-ordered. 
A necessary but not sufficient condition for stability 
is that all the main diagonal elements of the operations 
table be negative. 

Corollary 3.—The equilibrium is neutral if, and only 
if, the operations table forms a singular negative 
definite matrix. The determinant made up by its 
elements therefore vanishes. This is the essential 
criterion for the Determinant Method of evaluating 
buckling loads. It may be seen also from physical 
considerations that singularity of the matrix is a 
necessary condition for neutral equilibrium. In fact, 
when the structure buckles there must exist more than 
one deflection pattern which corresponds to equilib- 
rium—that is, the simultaneous linear equations 
represented by the operations table must have more 
than one solution. The equations are then dependent, 
the rank of the matrix of their coefficients is smaller 
than the order, the determinant of the coefficients 
must vanish, and the matrix is singular. 

Singularity of the matrix is a sufficient, as well as 
necessary, condition of neutral equilibrium only if it is 
stipulated that the matrix be negative definite. The 
next theorem, however, shows that in an investigation 
of the first buckling load, which is the only one that is 
of practical interest, there is no need of checking the 
fact that the matrix is negative definite. 


Theorem 2 


Let the loads P; throughout the structure be con- 
nected with a single variable load P by a known re- 
lationship, usually P; = K,P, where K, is a constant 
for any 7. Then, if the load P is increased from zero, 
the value of some principal minors of the matrix A 
will change; if the first minor to become zero is A,, 
then at the same time all minors with k or more rows 
will vanish. The load at which A; becomes zero is the 
first buckling load. 


A, ~ 0 may be assumed without loss of generality. 


Then, if A, is the first minor that vanishes, assume 
k <r, where r is the rank, and rearrange the rows 
(and columns, so as to preserve symmetry), as de- 
scribed in Table 2, so that A,-1Ay+1 < 0. The ele- 


ments of A, and therefore the determinants that contain 
them, are, in general, continuous functions of the load 
P, so that no determinant will change sign without 
first becoming zero. A,—, and A,+; are now of opposite 
signs; but certainly a load P can be found low enough 
so that the system is stable, the matrix nonsingular 
negative definite, and therefore A,-; and A,+, of the 
same sign, (k — 1) and (k + 1) being either both even 
or both odd. Therefore, either A,—; or A;,+; must 
have been zero before A,, and, as this is contrary to 
hypothesis, it is impossible that k <r. Then accord- 
ing to Table 2: A,+1 = ... = Ay-1 = A, = 0. The 
minors with less than r rows preserve their sign; 
therefore the matrix is singular negative definite, and 
the state of equilibrium is neutral. As this is, by 
hypothesis, the lowest load at which A, vanishes, 
the load is the first buckling load. 

The practical importance of this theorem is that, in 
evaluating the buckling load of a structure (Determi- 
nant Method), it is only necessary to find the load at 
which A, vanishes and then to check that A, has not 
vanished at any lower load. It should be noted that 
this check may not always be considered necessary and 
that it may often be obtained by an inspection of the 
functions making up the influence coefficients. 


Theorem 3 


Convergence to a finite unique set of values for all 
the unknown displacements in a straightforward re- 
laxation procedure is a necessary and sufficient criterion 
for the stability of the structure considered. 

This theorem is the basis for the Convergence 
Method of determining buckling loads. It differs 
practically only in nomenclature from theorem 1 of 
reference 11. Proof is omitted here, since it would 
follow closely that given in that reference (see also 
references 10 and 12). The proof is based on the fact 
that the total potential of the structure (to the negative 
of which the operations table corresponds) must be a 
minimum if the system is to be in stable equilibrium, 
a condition possible only when the work quadratic form 
based on the operations table is negative definite. 

Proof of this theorem for some particular relaxation 
procedures may be found in corollary 5. Work related 
to the Convergence Method may also be found in ref- 
erence 13. Some special examples in which a relaxa- 
tion process may converge even though stable equili- 
brium is not possible are discussed in reference 11 and 
corollary 5. 

It should be remembered that in complex problems 
it may be difficult to find a converging relaxation 
procedure even under stable conditions, but clearly 
in such cases divergence of the process should not be 
taken to imply instability. 


Theorem 4 


Every square submatrix of an operations table, taken 
so that its main diagonal is formed only by elements 
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chosen from the main diagonal of the complete matrix, 
represents a possible physical structure that will be 
called a substructure. Every substructure corresponds 
in fact to the original structure in which some of the 
points are rigidly fixed in some directions and the 
necessary reactions are there applied in those directions. 


Theorem 5 


All substructures of a stable system are stable, since 
their principal minors satisfy the requirements for a 
negative definite nonsingular quadratic form as given 
in Table 2. 

Corollary 4.—The lowest buckling load of any 
structure is either left unchanged or raised, but never 
lowered, by rigidly fixing some of its points. 


Theorem 6 


Equilibrium corresponds to a stationary value of the 
total potential V [Eq. (14)]. 

A typical one of the conditions for a stationary 
value of the function V may be written as 


OV /dx, = +2 ass) 


Eq. (18) is identical with the ith equation [Eq. (3)], 
represented by the operations table except for a factor 
of (—1), and is therefore one of the equilibrium condi- 
tions. It follows that, if at any load P the total 


‘potential V were approximately known in a particular 


problem in terms of several parameters, the nearest 
approximation to equilibrium would ‘be obtained by 
setting to zero the derivatives of V with respect to 
each of the parameters. 

Corollary 5.—As a corollary to theorem 6 a step- 
by-step method may be developed for obtaining 
approximate solutions of systems of linear equations, 
when approximate values for some of the unknowns 
are available or may be guessed. 

Break up the matrix A of the system into four sub- 
matrices, preserving symmetry, submatrices B and E 
being square, as follows: 


B Cc 
such that B and C include all the equations corre- 
sponding to the known displacements, and D and E 
those corresponding to the unknown ones. Let the 
known and unknown displacements be referred to as 
the B and E displacements, respectively. To set to 
zero the derivatives of the total potential V with 
respect to the E displacements with given values of 
the B displacements is equivalent (see theorem 6) to 
solving the equations represented by £, using as the 
right-hand side the negative of load terms P; plus the 
negative of the sum of the B displacements multiplied 
by the coefficients contained in D. This operation is 
the first cycle in the procedure. The process may be 
repeated now in a second cycle at the outset of which 


the newly found E displacements ‘are assumed and 
substituted in C and the B displacements solved for 
similarly to above. As many cycles may be used ag 
are needed to give the desired degree of approximation, 

Let the displacements x be considered as the inde- 
pendent variables and the total potential V as the 
dependent variable in Eq. (14). If then this equation 
represents a paraboloid, either an absolute maximum 
or an absolute minimum of the function V will exist, 
and all sections that may be cut through that surface 
by holding some of the variables constant will have a 


‘maximum or a minimum, respectively. The above 


step-by-step procedure will now be shown to converge 
in both of these cases. 

The total potential V, in the kth cycle, say, is based 
on some assumed E displacements, while the B dis- 
placements are found, according to theorem 6, so as 
to give either the maximum or the minimum value of V 
consistent with that assumption. Let these dis- 
placements be known as the FE, and B, displacements, 
respectively. The potential V,-: in the (k — 1)th 
cycle is made up of the E, displacements, while the 
B,-1 displacements are in general different from the 
B, displacements. The potential V, must then be 
nearer than the potential V,—-; to the absolute maxi- 
mum or minimum of the function V. Since this is true 
for any cycle k, the value of V will converge to the 
absolute maximum or minimum, respectively. 

The absolute maximum case corresponds to a quad- 
ratic form W, of Eq. (13), which is positive definite 
and therefore is related to an unstable system. The 
absolute minimum case corresponds to a quadratic 
form W,, which is negative definite and therefore is 
related to a stable system. As the two cases may be 
easily distinguished by the sign of the maiti diagonal 
elements (all positive and all negative, respectively), 
this procedure is seen to fall under the scope of theorem 
3. 

There remains now to investigate the case in which 
the function V has neither an absolute maximum nor 
minimum, but a saddle point. Two cases may be 
considered, in the first of which the function V varies 
monotonically with each cycle of the procedure pre- 
viously discussed and in the second of which the value 
of the function V oscillates. The former case will lead 
to convergent procedure only in exceptional examples, 
which, as suggested in reference 11, may be detected 
by carrying out a different relaxation process on the 
same problem. In the latter case the same test may 
be applied, or the variation of function V in each cycle 
may be checked. Oscillation of the value of this 
function in this relaxation process is sufficient to imply 
instability. 

A similar procedure may be developed by breaking 


up the given matrix in more than four submatrices, 


each submatrix not having necessarily the same number 
of rows. The special case in which each submatrix 
has only one row and one column corresponds to a 
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relaxation procedure in the Southwell method in which 


only individual displacements of one point at the . 


time are undertaken so as to balance the residual force 


existing at that point. 


Theorem 7 


If the quadratic form W, [Eq. (13)], representing 
the total potential when the terms P; are omitted, is 
evaluated for values of the displacements x corre- 
sponding to the deflected shape at the first buckling 
load, then W, will vanish if the system is in neutral 
equilibrium and will be negative if the system is below 
the first buckling load. 

The deflected shape at buckling is to be understood 
as consisting only of the distortions occurring at the 
moment of neutral equilibrium, which are in actuality 
superimposed upon the displacements existing imme- 
diately before buckling. Considering these latter 
existing displacements as a “‘particular’’ solution to the 
linear equations formed by the operations table in 
which the right-hand side loads P; are included and 
considering the superimposed buckling displacements 
as a “complementary” solution to the homogeneous 
(right-hand side zero) set of equations, it may be 
observed that the sum of the particular and comple- 
mentary solutions constitute the total distortions 
immediately after the moment of buckling and again 
satisfy the equilibrium conditions as specified by the 
operations table. This clarifies the omission of the 
external loads for the calculation of the total work in 
the statement and proof of this theorem. 

Proof may be obtained directly from theorem 1. 
It may be noticed that at the moment of buckling all 
the equations exemplified by Eq. (3) will be satisfied 
by the buckling distortion pattern when the terms 
P;, are set equal to zero. Since W; is in reality the sum 
of the left-hand side terms in that equation, multiplied 
for each value of 7 by a constant, the value of W, will 
be zero. Below the first buckling load the structure 
is stable; hence, the quadratic form is nonsingular 


' negative definite, and its value will be negative for any 


set of displacements (such as the buckling displace- 
ments) which are not zero throughout the structure. 
This theorem will also hold if the quadratic form Q 
of Eq. (8) is substituted for W,. 

Theorems 6 and 7 are the basis for the Energy 
Method of calculation of buckling loads. 


THREE METHODS FOR THE NUMERICAL DETERMINATION 
OF BUCKLING LOADS 


Three methods for the determination of buckling 
loads are outlined in this section, proof of their validity 
being contained in the previous theoretical considera- 
tions. Illustrations of these procedures are given in 
the last section of this paper. 

Each method requires a knowledge of the stress 
distribution throughout the structure at buckling. 


This may be most conveniently obtained by setting 
up and solving the operations table, based on an 
assumed applied load. In fact, some of the terms 
appearing in this operations table will be used in the 
next step of each method. This next step is the 
setting up of the operations table for the problem of 
instability, which differs from the previous operations 
table inasmuch as in the influence coefficients terms 
must be included which are functions of the applied 
load P. 

Determination of the critical value of P may then 
be carried out by one of the following methods: 


(1) Determinant Method 


Evaluate the determinant formed by the elements 
of the operations table for the buckling problem; 
set the resulting expression equal to zero; calculate 
the load P from the equation obtained. The lowest 
root of this equation is the first buckling load. In 
many cases, however, the operations table may involve 
too complicated trigonometric functions to allow such 
a direct evaluation of the unknown load P. An 
alternative procedure may then be followed, which 
consists in assuming several values of P and calculating 
the value of the determinant corresponding to each 


(which is now a purely numerical process). A guide . 


in the assumption of the loads P is given by Corollary 2 
in the section on Theoretical Considerations. If the 
determinant value is plotted against the loads P, the 
buckling load P,, can be read off from oné of the 
intercepts of the resulting curve. It might be con- 
venient at times to reduce the order of the determinant 
by operations on rows and columns before applying 
the above process. Evaluation of determinants is 
simplified at times by reduction to a triangular form 
by a method such as that of reference 7. 


(2) Convergence Method 


Any convenient straightforward relaxation pro- 
cedure may be employed, the following one being 
perhaps the most obvious. Set up the operations 
table for a given value of the load P. Assume a unit 
load to be acting at one of the points in the structure 
which is likely to be displaced during buckling. Bal- 
ance this load by a motion of this same point, thus 
causing residual forces to exist at the others. These 
forces may now be relaxed by motions of all the points 
except the original one, until the latter is the only one 
on which a residual is left. Then, according to the 
convergence criterion, if this new residual is smaller 
than the original 1-lb. load, the system is stable— 
that is, the load P was assumed too low. If it is larger 
than 1, it is unstable—that is, the load P was assumed 
too high. If it is equal to 1, it is neutral—that is, 
assumed load P equals P.,.. 

The relaxation procedure should be carried out for 
several values of the applied load P, with the residual 
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obtained in each case plotted against the load P. The 


intercept of the resulting curve on the line of unit. 


residual is the critical load P.,.. 


(3) Energy Method 
Assume values for the displacements of some of the 


‘points in the structure at buckling. An often con- 


venient guide to this assumption may be given by the 
deflected shape at buckling obtained from the solution 
of the same problem, in which, however, some points 
or members have been neglected or some similar 
simplifying assumption has been made. Calculate 
the quadratic form Q [Eqs. (7) or (8)] for an assumed 
value of the load P on the basis of the assumed de- 
flections, differentiate it with respect to the displace- 
ments that still remain unknown. Setting these 
differential coefficients equal to zero leads to a system 
of linear equations from which the remaining unknowns 
may be found. The numerical example of the last 
section shows how these equations may be obtained 
directly from the operations table. Substitute into 
the expression for Q the values found for the unknowns 
and obtain a value for Q. A negative value of Q 
indicates that the assumed load P is too low; a positive 
value, that it is too high. Repetition of the entire 
procedure for other values of the load P yields a curve 
representing the variation of Q with P. The value 
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of P at which Q = 0 may then be read off and is the 
critical buckling load P,,.. 

It may be noted that all three methods require the 
assumption of various loads P. A guide in this as. 
sumption may be often given by a knowledge (see 
Table 2) of the signs of the determinant formed by the 
elements of the operations table, of the main diagonal 
elements, and of all the other principal minors. For 
example, should any main diagonal elements in an 
operations table be zero or positive or any two-rowed 
principal minor zero or negative, or should the complete 
determinant be positive if the number of rows is odd 
or negative if the number of rows is even, it will imme- 
diately be concluded that the system is unstable and, 
consequently, the assumed load too high. 

It should be remarked that, if the buckling load 
obtained differs materially from that originally assumed 
in the calculation of the stress distribution, the entire 
procedure must be repeated on the basis of the buckling 
load obtained. 

In order to ensure that the first and not a higher 
buckling load was found from the calculation, it may be 
necessary to check the stability at loads lower than the 
buckling load found. 


EXPERIMENTAL INVESTIGATION 


Two identical specimens of sheet, stringer, and ring 
combinations were tested in PIBAL, so as to obtain 
an experimental check on the above-described methods 
of determining buckling loads. The same structure 
is used in the numerical examples of the next section. 

A drawing of the structure is shown in Fig. 3; photo- 
graphs of a specimen after buckling are presented in 
Figs. 4 and 5. As the photographs indicate, the ends 
of the rings and sheet panels were bolted to a rigid 
welded steel framework. This structure and _ the 
attachment of the specimens to it were considered 
sufficiently rigid so as not to distort under load. 

Each end of the stringer extended beyond the end 


rings and was clamped between two hardened steel 


grip fittings serrated in order to prevent slippage and 
squared up parallel to the end surfaces of the testing 
machine. This type of attachment was expected to 
give a rigid end fixation for the stringer. To check 
the validity of this assumption, a test was run on a 
straight 24S-T aluminum-allov har of */s-in. square 
section. Two of the fittings were attached to each end 
of the bar, and an axial compressive load was applied 
through them in the same manner as to the sheet and 
stringer combinations. The experimental failing load 
checked within less than 1 per cent with the Euler load 
calculated according to fixed end conditions. This 
was considered sufficient proof of the rigidity of the 
fitting attachment. 

The load was applied to the sheet and stringer 
combination specimens by a 200,000-lb. capacity 
Riehle Brothers lever-type testing machine. The 
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Fic. 4. Back view of test specimen 2 after buckling. 


strain in the stringer of each specimen was measured 
by pairs of Baldwin-Southwark SR-4 type A-1 wire 
resistance strain gages placed in the middle of each 
field. The strain distribution obtained is discussed 
in the next section together with the comparison with 
theoretical results. 

The applied loads P at buckling for the first and 
second specimens were 3,820 and 3,780 ibs., respec- 
tively. The deflected shape of the specimens at 
failure may best be seen from Fig. 5. The stringer 
buckled inwardly in the end fields and outwardly in 
the middle fields. The rings were also distorted from 
their original circular shape. 

Good agreement between results obtained with the 
two specimens was found with regard to strain dis- 
tribution, buckling load, and deflected shape. 


NUMERICAL EXAMPLES 


The buckling load for the structure of Fig. 3 is 
calculated in this section by means of each of the pre- 
viously described methods. It is first necessary to 
determine the stress distribution in the structure, 
and this may be done by means of an operations table. 

The operations table is set up as described in ref- 
erence 6, and takes into account the symmetry of both 
loading and structure. Point A, in reality, could not 
move in the radial direction, while the rigid steel frame- 
work attached to the ends of the rings was actually 


eal 


Fic. 5. Side view of test specimen 2 after buckling. 


free to move as a rigid body. As in the operations 
table, point A is considered movable and the ends of 
the rings are considered fixed, the magnitude of the 
displacements obtained from the solution of the 
operations table should be interpreted as relative to 
the plane of the ends of the rings. The operations 
table follows as based on a unit displacement of 
0.001 in.: 


TA xB XA —R.HLS. 
Rg —3.67844 0.01184 0.3704 0.3704 O 
Ra 0.01184 —3.67844 — 0.3704 —0.3704 
Xp 0.3704 —0.3704 — 520.8 148.6 0 
Xa 0.3704 —0.3704 148.6 -—171.8 P 


In this operations table r, R, x, and X represent the 
radial displacement, radial force, axial displacement, 


and axial force, respectively, at the point indicated by 


the subscript. 

The influence coefficients depend upon the value of 
the effective width of sheet acting with the stringer 
and the value of the effective shear modulus Gy. 
for the sheet panels. The following values were used 
for the area A.q and the bending rigidity EJ of the 
stringer plus the effective width of sheet: 


fields AB or A’B’ Ag, = 0.1476 sq.in. 
EI = 19,743 Ibs. in.? 
Ges./ Go = 0.202 (20) 
field BB’ = 0.16 sq.in. 
EI = 22,032 Ibs. in.? 
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The effective area was chosen smaller in the end panels 
than in the middle section, because in the former the 
strain near the stringer was expected to be higher and 
because it was thought that the sheet was not likely 
to be fully carrying in that region since in the experi- 
ments the load was applied only to the stringer. No 
values for the effective shear modulus of curved panels 
buckled because of shear and compression could be 
found in the literature, but, according to data obtained 
from flat panels and circular cylinders in torsion™ 
and combined torsion and compression,* the assumed 
value does not appear unreasonable. 

Solution of the operations table by the method of 
reference 7 leads to the following values for the un- 
known displacements, in inches: 


0.0077353P X 10-* 


x = 

Xp = 0.0022085P X 10-* (21) 
ra = —0.0009981P X 10-% 

tp = 0.0009981P xX 10-* 


The force P’ in each bar may now be found from 
the formula 


1947 


= (AL) EA (22) 


where AL is the shortening of the bar and is equal to the 
difference of the displacements of the two ends of the 
bar. In this case 


Pap = 0. 
Psa = 0.76978Pf 


while the corresponding experimental values are 


average of Py, and =.0.873P 
Pop = 0.765P ( ) 


(23) 


The theoretical force distribution of Eq. (23) may 
now be used in the calculation of the buckling load. 
Using again the values given in Eq. (20) for effective 
areas and. shear modulus, neglecting the axial dis- 
placements (which were found in several other examples 
not to influence appreciably the critical load), consider- 
ing the bending of the stringers (Table 1) and torsior. 
of the rings [Eq. (6)], and assuming symmetry of 
loading, structure, and distorted shape, one << 
establish the following operations table: 


Ra Re Me 
YA —3.67844 — (Pks/D,)1 0.01184 + (Pks/D;)1 — 
0.01184 + (Pks/D;)1 —3.67844 — (Pks/D,)1 —c)/Di)1 
wR — ©)/Di)1 (PQ — ¢)/Dij1 0.4 — [P(s — kLc)/(kDi)]1 + [P(RL — 


— kLe) (EDs) 11 


The subscript J for an influence coefficient refers to 
the end fields, the subscript JJ refers to the middle field, 
and the load P is given in units of 1,000 Ibs. The 
symbols wz, and M, refer to rotations and moments, 
respectively, at point B in the direction indicated 
vectorially in Fig. 3. 


—1.84514 — (Pks/D, [Pa 


The value of | A| was calculated from Eq. (25) for four 
different values of the load P, and the results were 
plotted as shown in Fig. 6. The intercept at which the 
determinant vanishes is the critical load and was 
read off as 


P.,, = 4,008 Ibs. (26) 


‘ 


(2) Convergence Method 


Relaxation of a 1-lb. force acting at point B in the 
radial direction was carried out for four different 
values of the load P on the basis of the operations 


From this operations table the critical load P.,, will 
be now obtained by each of the three previously out- 
lined methods. 

(1) Determinant Method 

By operations on rows and columns the determinant 
|A| represented by the operations table may be 
reduced to 


— [P(2s — kL — kLc)/(kD,))11| (25) 


table given above. It is shown here in Table 3 for 
one such load chosen in the stable range. 

The last two displacements undertaken in the relax- 
ation section of Table 3 were obtained by multiplying 
the two preceding motions by a factor of (0.21605)+ 
(0.49086 — 0.21605). 

When the relaxation was repeated for a load P = 
4,060 Ibs., the new residual obtained was 29.4465, and 
it was therefore concluded that at this load the struc- 
ture is unstable. A rapid diverging procedure is 
therefore encountered for a load only about 50 Ibs. 
higher than the critical load. 
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The new residual may be plotted as in Fig. 7 against 
the corresponding load, and the intercept on the line 
of new residual = 1 may be read off to be 


P.., = 4,008 Ibs. (27) 


(3) Energy Method 


An approximation to the distorted shape of the struc- 
ture at buckling will be obtained from the solution 
of the problem which has been simplified by neglecting 
the bending of the stringer. With the aid of Table 1, 
if P is expressed in units of 1,000 Ibs., the operations 
table for this case is found to be: 


TA 
Ra —3.67844 + P/(9.64) 0.01184 — P/(9.64) 
Rp 0.01184 — P/(9.64) —3.67844 + P/(9.64) 


By the determinant method, the value of the load P 
which makes this two-rowed determinant vanish may 
be found from direct evaluation to be P = 17,790 Ibs. 
Substituting this load back into the influence coeffi- 
cients, the linear equations represented by the opera- 
tions table become , 


6 
72 
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APPLIED LOAD P, LB. 
Fic. 6. Numerical example for the determinant method. 


TABLE 3 
P = 4,000 Ibs. 
Operations 
Forces ~ 
Displacements Ra Me 
—3.48956 —0.17704 — 0.86725 
—0.177 —3.48956 0.86725 
WB — 0.86725 0.86725 — 0.48972 
Relaxation 
Ra Re Ms 
Assumed Forces 0 1 0 
rg = 0.28657 — 0.05073 —1 0.24853 
wg = 0.50749 —0.44013 0.44013 — 0.24853 
z= —0.49086 0.44013 0 
ra = —0.14067 0.49086 0.02490 0.12200 
wg = 0.24912 — 0.21605 0.21605 —0.12200 
—0.21605 0.61809 0 
ra = —0.11059 0.38591 0.11048 0.09591 
wp = 0.19585 —0.16986 0.16986 —0.09591 
New residual = 0 0.89843 0 


0.89843 < 1 therefore converging, therefore stable. 
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Fic. 7. Numerical example for convergence method. 


4050 


—1.8333r, — 1.8333r, = 0 \ (28) 


—1.8333r4 — 1.8333r, 0 


Il 


The solution to these equations is the approximate 
buckling distorted shape and is 


= (29) 

The energy function Q corresponding to a load P = 
4,000 Ibs. (see the operations section of Table 3) is 

Q = —3.48956r4? + 2(—0.17704)rare + 
2(—0.86725)rawg — 3.48956rz,? + 

2(0.86725)rgwg — 0.48972w,? 

With the assumption of Eq. (29) this reduces to 

= —6.625 + 3.469(wp/rp) 0.4897 (wp/rp)?* (31) 

Setting equal to zero the derivative of this expression 

with respect to (wg/rg) one obtains 

(wp/rp) = 3.542 

This value might have been found, perhaps more 


easily, by considering the equation in the operations 
table which corresponds to w,—namely, 


—0.86725r, + 0.86725r_g — 0.48972w, = 0 


(30) 


(32) 


(33) 

If the assumption of Eq. (29) is used, the value of 
(wp/rg) from Eq. (33) is: 

(—0.86725)(—1) + 0.86725 
= 
0.48972 
Substitution of the value of wz/rz in Eq. (31) yields 
Q = —0.4815r,? (35) 


A negative value of Q indicates that the system is 
stable under a load of 4,000 Ibs. 

Evaluation of the quadratic form Q was carried out 
for four different values of the load P; the value of 


3.542 (34) 
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of the 
28 
(23) | 
| 

may 
ective 
sider- 
rsior. 
ry of f 
may 
will 
out- 
lant 

for 
) + 
nd 

is 
DS. 


348 JOURNAL OF THE AERONAUTICAL SCIENCES—JUNE, 


6 

4 

‘a 

=4008 LB, 


APPLIED LOAD P, LB. 
Fic. 8. Numerical example for the energy method. 


Q was plotted against the load P in Fig. 8; and the 
Q = 0 intercept was read off as 


Por, = 4,008 Ibs. (36) 


Of course, the accuracy of the Energy Method 
depends on the original assumption for the deflected 
shape. Any variation in the assumed shape from the 
true buckling distortions will result in a value for the 
buckling load higher than reality. The method, 
however, is not believed to be too sensitive to changes 
in deflected shape. In fact, let the assumption of 
Eq. (29), which happens to coincide with the true 
distortions, be replaced by 


= —0.75rg (37) 


that is, by an arbitrary error of 25°per cent; the 
buckling load then obtained is 


P.,, = 4,014 Ibs. (38) 


which is practically the same as that previously 
calculated. 

It may be difficult in some problems to know whether 
a good approximation to the buckling load has been 
obtained from this method when a deflected shape has 
been assumed at the outset. The accuracy may be 
improved by decreasing the number of points of which 
the deflection is assumed. A buckling load closer to 
reality will also be obtained if the calculations are 
repeated on the basis of assumed deflections that co- 
incide with those previously found by minimization of 
the quadratic form Q and if the deflections considered 
as unknowns are those occurring at the points where 
the deflections were assumed at the start. 

The results obtained with the Determinant and Con- 
vergence Methods are, as they should be, identical 
with those obtained from the Energy Method in which 
the true buckling distortion pattern [Eq. (29)] is used. 
The theoretical buckling load differs from the empirical 
value by about 5.5 per cent of the experimental average. 
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Discussion 


Mrs. Hilda Geiringer, Professor of Mathematics, Wheaton 
College: ‘This paper presents three efficient numerical methods 
for the solution of buckling problems. The theoretical back- 
ground of the computations is formulated in a series of interesting 
and useful theorems of a mechanical nature. It is attempted 
to comment briefly on some of them from a mathematical point 
of view. The author’s investigation is based on Southwell’s 
relaxation method which has been successfully applied to various 
problems in mechanics." 

“Consider the example of a redundant framework often used 
by Southwell. The configuration of the elastic structure is 
specified by n generalized coordinates x, x2, and x, which de- 
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CALCULATION OF ELASTIC INSTABILITY 


termine the displacements of the joints from the position of an 
yncharged equilibrium. The Maxwell equations of the system 


are of the form 


+r=0, = ani (1) 


where the r; are external forces and the coefficients a;, are influence 
coefficients in a generalized Hooke’s law, which are symmetric, 
according to Maxwell’s reciprocal relations. Following Dr. 


Boley, denote by W; = (?/ 2) Doauntcte the internal work and by 
i, 


= —W, the strain energy of the structure, which is positive 
definite if the position of equilibrium is stable. The total po- 
tential V—i.e., the sum of the potential of the external forces 


and the total strain energy is then 


V = —([('/2) + Pixi), Gin-= ani 


Avoiding the negative sign, one may put 


V= + 7iXi, ik = am (2) 
i, 
The ith equation [Eq. (1) ] follows from Eq. (2) by differentiating 
V with respect to xj. 

“It has been pointed out, occasionally, that the relaxation 
method for the solution of Eq. (1) is based, mathematically, on 
the iteration method, known by the name of Seydel* which may 
be described as follows. Assume a first set of arbitrary approxi- 
mation values x, ...x,“, and correct the first value by solving 
the first equation with respect to the first unknown, thus deter- 
mining a new value x;; next a value x“) is found from the 
second equation using x, x, ...x,{, ete. Finally, 
follows from the nth equation by means of x, x3, ...%n-1. 
Anew cycle starts where x; is computed from the ith equation 
by means of .x;-1, xi+1, ...2%n, ete. It has been 
shown in 19294 that this procedure converges if the matrix A of 
the ai, is positive definite. This is the well-known condition for 
stable equilibrium (see the author’s Theorem 1). 

“The described iteration procedure constitutes a relaxation 
procedure, since the ‘successive steps in this iteration may be 
interpreted as a systematic relaxation of the constraints repre- 
sented by the residual forces.’ In fact, if the vector xi™ is 
different from the solution x; of Eq. (1), the left side of the 
ith equation [Eq. (1)] is not equal to zero but to a ‘residual 
force’ R;\”): 


= p> ai + % (3) 
Jj 


and the iteration consists in correcting the x; so as to reduce 
the residual forces to zero. [See reference 5 and: ‘We whittle 
away at the given data until every residual quantity has been 
brought within a certain margin of uncertainty’ (Southwell?) ]. 
It is important for the practical application that—as can be 
concluded from the convergence proof—the order in which the 
Eqs. (1) are used in the above iteration is irrelevant for the con- 
vergence. As long as each equation is used again and again, 
they can be used in any order. The iteration may, e.g., be 
worked so as to correct each time that unknown x; for which 
the corresponding residuum R;™ is the greatest (Southwell). 
Another noteworthy fact is that, in case of a positive definite 
matrix A, the starting point x;“) may be chosen quite arbitrarily. 

“Theorem 3 of Boley (stated before by Professor Hoff* 7) 
implies, if I understand it right, that positive definiteness of A 
is not only sufficient but also necessary for the convergence of the 
iteration. Trying to prove this, one may reason as follows. 
Denote by AV the difference of the V-values [defined in Eq. 
(2)] before and after a single step in the iteration: 


AV = V(x Ot), 


xa) (4) 


xt), 


349 


It can be shown easily that 
AV = — = (4’) 


where R,) has been defined in Eq. (3). [Eq. (4’) constitutes 
the essential point in the proof of the sufficiency of positive 
definiteness for the iteration.‘] Now assume that the elements 
aiz Of a nonsingular matrix A are symmetric and ai > 0. If A 
were indefinite, there would exist an initial value xi) such 
that the corresponding 3 = V(x, ... x,) is smaller than 
zo = Vix ... xn), where x; denotes the solution of Eq. (1). 
But since, because of Eq. (4’) and ay > 0, each single step de- 
creases the value of V, the point 2 could never be reached. 
(There may, of course, still be convergence for a particular set 
of starting points.) We thus formulate what might be essen- 
tially the algebraic equivalent of Theorem 3: 

“(A) Let OV/Ox; = 0 with V defined in Eq. (2) be the ith of 
the Eqs. (1) and use in the above explained iteration method the 
ith equation to correct the value of the ith unknown (i = 1,...n). 

“Then: 1) If Q = Dainxixe is positive definite the iteration 
converges for any starting point and for any order of the equations 
2). If (a) the determinant A # 0, (b) aik = ani, (c) ax > O, 
and if there is one order of equations such that the iteration converges 
for that order and for an arbitrary starting point, the quadratic 
form is positive definite. (Consequently, the procedure then 
converges for any order and for any starting point.) Let us note 
however, that the procedure may converge also for indefinite 
Q for certain starting points. 

“It may be of some interest in this connection to state, for 
general az, the necessary and sufficient conditions of convergence 
of the iteration in question. We may assume that the equations 
are such that the unknown x; really occurs in the ith equation, 
i.e., that aii # 0; we may then assume without loss of generality 
that aij = +1. The iteration procedure is defined by: 


+75 = 0 (5) 
If we introduce the error vector 
2” = x” — x, (6) 


where x; is the solution of Eq. (1), we get from Eq. (5): 


+ as zo”) +... + aintn”) = 0 
+ +... + == 0 
+... + = 0 


Then the following theorem holds which is a consequence of 
well-known algebraic theories. 

“(B) The iteration method (relaxation method) defined in 
Eq. (5) converges toward the solution of (1), > O0asy 
if and only if all roots pi of the equation 


Qin 
anp p cee Qin 


G.(p) = = pF,-1(p) = 0 (7) 
... p 
are less than one in absolute value. 
“If, however, px = 1 is an r-tuple root of F,-:(p) = 0 and all 
other roots are less than one in absolute value the iteration (5) with 
1% = O converges toward a solution of the homogeneous problem 


Drains = 0 if and only if the rank of A equals n — r. 


“It can be seen that the error vector 2;”) tends toward zero 
for every choice of starting point x; if all roots p; of Eq. (7) 
are less than one in an absolute value. If, however, some of these 
roots have absolute values = 1, there may still be convergence for 
a certain set of starting points xi, This is worth mentioning 
from the point of view of the practical application: If it is desired, 
as in Boley’s approach, to use convergence as a criterion for 
positive definiteness, i.e., for stability it must be remembered 
that there may be convergence of the iteration for infinitely many 


Theory, 
col Pre 
es. 
Strip 

| 
946. 
ton 
ods 
ck- ie 
ing 
ted 
int 
ll’s 
ed 
is 

. 


350 


starting points in spite of an indefinite matrix A. This possibility 
is indicated in the author’s paper. On the other hand in case 
of a positive definite matrix A, the roots of Eq. (7), and of the 
(n — 1)! similar equations that correspond to another order of 
Eqs. (1), are all less than one in value. 

“A further remark may be made regarding Corollary V of 
the author’s paper. For the application of the relaxation method 
it need not be assumed that each equation be used separately. 
‘Blocks’ of equations may be used to correct simultaneously 
corresponding groups of unknowns. This seems to correspond 
algebraically to ‘group iteration.’* The convergence of this 
procedure was established in reference 4 in case of a definite 
matrix A. In the present paper the author considers group 
iteration for the indefinite case as well, stating as a result that 
oscillation of the value of V under this procedure is sufficient for 
instability of the system. 

“In ‘group iteration’ the ‘single step’ consists in solving 
exactly a small group of equations with respect to a corresponding 
group of unknowns. We may, e.g., for m = 8 split the equations 
into three groups, the first consisting of the first two equations 
with the first two unknowns, the second consisting of the following 
three equations with the corresponding unknowns in the main 
diagonal, the last group comprising the last three equations with 
the remaining unknowns. It must be assumed that the re- 
spective coefficient determinants of second, third, and third 
order are not zero. Starting with a first set x;“(i = 3,...8), 
an approximation x, x2 is found by solving simultaneously 
the first two equations. Then, using x, x2 and x6, x7, 
xs, we find x, x;@ by solving simultaneously the 
following three equations, etc. Then the same scheme starts 
anew, etc. 

“The necessary and sufficient condition of convergence 
for every starting point can be written immediately. Consider, 
e.g., 2 = 5 equations with groups (1, 2), (3, 4), and (5). Then 
Eq. (7) is to be replaced by the equation 


=| ano ant ano ax | = 0 (8) 


“(C) The necessary and sufficient condition of convergence, for any 
starting point, is that all roots of Kn-2(¢) = O are less than one in 
absolute value with an additional statement (analogous to the previ- 
ous one) in case of roots equal to one. Again, for certain sets of 
starting points convergence may happen if some of the Eigen 
values are 21. 
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“Group iteration may prove advantageous if certain groups of 
variables associated to certain groups of equations form a natural 
mechanical unit and if perhaps the solution for such a ‘block’ js 
readily available; or it may be combined with the use of a com- 
puting machine for finding the exact solutions of the small groups 
of equations. (See also a comment by Geiringer, Vol. 48, Bull. 
Am. Math. Soc., p. 370.) 

“If A is positive definite, the determinants of the various groups 
of variables, which are principal minors, are certainly not zero 
and group iteration converges for any starting point. Moreover, 
it follows easily that convergence will be accelerated by group 
iteration. More precisely: 

“(D) In case of @ positive definite matrix A, a ‘step’ in group 
iteration applied to a group of s > 1 equations and unknowns re- 
duces the value of the potential V in Eq. (2) by more than s suc- 
cessive single steps applied to the same group of equations would do, 
Geometrically, this is fairly obvious and follows easily by com- 
putation as well. 

“In the mechanical ‘application to an n-fold indeterminate 
structure each of the single steps in a single-step iteration means 
that we solve a statically simply indeterminate system, thus re- 
ducing the solution of an n-fold indeterminate system to the re- 
peated.solution of simply indeterminate systems.* In this sense 
group iteration means that instead of reducing the solution of 
the n-fold indeterminate system to the repeated solution of simply 
indeterminate systems we use as intermediate steps s-fold 
indeterminate systems. The exact solution of such an s-fold 
indeterminate subsystem with respect to the corresponding group 
of unknowns may be worked out in amy convenient way, even 
by means of equations that are not the original Maxwell equa- 
tions. 
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SUMMARY 


The pressure distribution over an oscillating, thin airfoil of 
infinite span advancing at a supersonic speed is calculated by 
integrating the disturbances due to an appropriate distribution 
of sources over the forepart of the Mach cone for the point at 
which the pressure is desired. The result is expressed in terms 
of a Green’s function giving the pressure at any point on the 
airfoil due to unit harmonic velocity at any other point on the 
airfoil. The lift and mid-chord moment due to plunging and 
pitching are expressed in terms of four complex coefficients 
Aga, Aah, Aca, and A¢p which coefficients are functions of the 
reduced frequency (&) and the Mach Number (1). Asymptotic 
results are given for both small k and large M@. The results are 
plotted for M = 1.2, 1.4, 1.6, and 2 and fork = 0 to 1.5. They 
are found to agree with the previous work of Possio’ (nu- 
merically) and Borbely® (analytically). The results of Theodor- 
sen! for M =0 and those of Possio* for 1 = 0.7 are given for 
purposes of comparison. 


NOTATION 


amplitude of plunging or pitching motions 


Thy a 

b = semichord 

c = velocity of sound 

h = plunging motion 

j=(-1)”? = imaginary-unit 

k =(wb/U) = reduced frequency 

p = pressure 

q = vector velocity 

r = radius 

t = time 

w = vertical velocity at point on airfoil 

& 9% = coordinates (wind along x, airfoil span along 
y, 2 positive up) 

Ay = dimensionless flutter coefficient giving force 
in ith degree of freedom due to motion in 

: jth degree of freedom 
G(x, ¢) = Green’s function giving potential at.x due to 
velocity at 

Jo = Bessel function, first king, zero order 

= lift 

M = Mach Number (l//c), also mid-chord. stalling 
moment 

U = free-stream velocity 

a = pitching motion 


Received November 5, 1946. 

* Since submitting the manuscript for the present paper, the 
author has obtained a copy of N.A.C.A.T.N. No. 1158, “Flutter 
and Oscillating Air-Force Calculations for an Airfoil in a Two- 
Dimensional Supersonic Flow,” by I. E. Garrick and S. I. Rubi- 
now, issued in October, 1946. This report gives tabular results 
that are equivalent to those of the present paper for M = 10/9, 
5/4, 2, 5/2, 10/7, 5/3, 10/3, and 5 and for values of K ranging 
from 0.01 to 10. The application of these results to certain 
cases of bending-torsion flutter is also considered. 
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The Aerodynamic Forces on an Oscillating 
Airfoil at Supersonic Speeds’ 


JOHN W. MILES+ 
University of California at Los Angeles 


(M2 — 1)'72 


8 = 

tn = source coordinates on (x, y) axes 
kK = kM(M? — 

= (1 — 

= velocity potential 

6 = 

@ => 


angular frequency 


INTRODUCTION 


— PROBLEM TO BE SOLVED consists essentially 
in finding the Green’s function G(x, &), giving 
the pressure p(x, ¢) at the point x on the airfoil shown 
in Fig. 1, due to unit normal harmonic velocity e* 
at the point ¢. The airfoil occupies the region bounded 
by x = ¥1, y = + © and is considered to be sufficiently 
thin and to produce sufficiently small perturbation 
velocities to allow the specification of the normal 
velocity in the plane z = 0. The free-stream velocity 
U = Mc is in the direction of the positive x axis. 

In what follows, } will be taken as the unit of length, 
where 2b is the chord of the airfoil. 


EQUATION OF MOTION 


Assuming that the perturbation pressure due to the 
airfoil oscillation is small compared to the free-stream 
pressure, it follows from linearized acoustical theory? 
that the vector velocity g may be specified by the 
gradient of a scalar potential ¢: 


G(x, 2, 1) = Vo(x, 2 (1) 


It also follows that ¢ satisfies the Helmholtz equation, 
remembering that differentiation with respect to time 
must be effected by the operator 


)/d] + )} 


X=#/ 


Fic. 1. Section of two-dimensional airfoil illustrating positive 
deplections in plunge (h), pitch (a), and elevon flapping (8). 
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since it is desired to describe the motion of a point 
fixed on the airfoil (and not fixed in the fluid). Linear- 
izing this operator by assuming that the velocities due 
to the airfoil are small compared to U then gives 


A more detailed derivation of Eq. (2) is given in ref- 
erence 5. 

In order to solve Eq. (2) it is expedient to introduce 
the transformation* 


ux (3a) 
= (1/u)t + w(U/c?)x (3b) 
= (1 — (3c) 


which accomplishes the desired result of transforming 
Eq. (2) to 


[V2,, — = 0 (4) 


and yet leaves (0/0/) = y(d0/d2), so that the ad- 
vantages of harmonic time dependence are retained. 


SOLUTION FOR A POINT SOURCE 


A fundamental solution to Eq. (4) representing a 
source at 7= 0 is 


= (5a) 


r= + y? + 2? (5b) 
where f[? * (7/c)] is an arbitrary function. Since it 
is desired to have ¢ become infinite as r~! at r = 0 
to represent an outward traveling disturbance vanishing 
at ry = o and to exhibit harmonic time variation, the 
function f is chosen as 


fl — = (6) 


Transforming back fo the (x, y, 2, 4) coordinates and 
dropping the e’ factor, it is seen that the potential 
due to a source of strength —q(£, )didn at the point 
(é, 7, 0) is given by 


u = [(x — — — n)? — (6b) 
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By differentiating Eq. (6a) with respect to 2, it is 
found that w = 0¢/0z vanishes everywhere in the 
plane z = 0 except at the point (£, , 0)—i.e., at the 
source—where it is indeterminant. To evaluate 
$(x, y, 2) in the vicinity of (&, », 0), it is necessary to 
divide the source area didn into smaller elements 
dét’dn’ and to integrate the contribution of these sub- 
divisions; moreover, since the variation of the ex- 
ponential is negligible compared to the variation of 
1/u over the region dédn at (£, n, 0), it is permissible 
to take the exponential outside of the integral sign. 
Hence, the local problem is the same as the non- 


oscillatory case (w = 0), and the integration for this 


case has been carried out by Puckett’ with the result 
= (1/m)w(E, 0) (8) 


SOLUTION FOR A LINE SOURCE 


In the two-dimensional problem it is assumed that 
w is independent of y and may, therefore, be written 
w(t) in the z = 0 plane. It is convenient to choose 
the representative section of the infinite span airfoil 
in the plane y = 0. 

Since the potential at any point is determined only 
by those sources in the fore part of the Mach cone 
extended from that point, the potential at the point x 
on the representative section is determined only by the 
velocities of those points lying in the area on the wing 
ahead of x which is cut off by the Mach cone—ice,, 
the (, area bounded by = — &)/8] and 
x = —1, where B = (M? — 1)". Hence, the potential 
at x is given by 


G(x, §) = (10) 


To evaluate Eq. (10), the transformation sin 9 = 
[8n/(x — £&)] is introduced; substituting u from 
Eq. (6b), where y = z = 0, yields* r 


Mk wh 
= ——,k=—, = M*—1 7 
whence 
* Eqs. (3) represent the superposition of a Lorentz transforma- ‘ad 
tion and a transformation bringing the airfoil to rest in the flow. G(x, &) = e”* oe Jo[x(x — &)] (12) 
PRESSURE ON THE AIRFOIL 


Bernoulli’s theorem for the present problem may be written as* 


(Op/dt) + 1/2(V4)? + Sdp/p = fit) (13) 


Th 
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Linearizing, the pressure is given by 

p(x) = —p[U(0g/dx) + (0¢/0F)] (14) 
Carrying out the required differentiations, remembering the (0/0f) = jw for the harmonic motion considered, 
the result is 


x 0G: . x 
The second term in Eq. (15) may be integrated by parts to yield 
U x ow x 


The result, Eq. (16), is valid (assuming linearized flow) for any two-dimensional thin airfoil in the supersonic 
régime. In order to proceed further, it will be assumed that the airfoil is a flat plate and that the motion is either 
plunging (h) or pitching (a) about the mid-chord* as defined by the equations 


z,(x, t) = (17) 
t) = —a,bxei (18) 
The velocities due to these displacements are given by 
w(x, t) = + U(dz/dx) (19) 
Substituting in Eq. (16) yields ; 
‘jup 1+x 
p,(x, t) = ( + 1) + jk (20) 
j 1+<x l+<x 1+x 
Pal(x, t) = + +1) -—2 G(v)dv — jkx i G(v)dv + jk Good fom (21) 
us 


The lift and moment about the mid-chord due to the 
pressure p,(x, ¢) are given by 


L,(t) 


= =2 t)dx (23) Ay = Ry + Glin Se = Ra + (29) 


where a positive moment acts to raise the leading edge. An = 4Ro + (2/k)Io — 2Ri (30) 
Substituting Eqs. (20) and (21) into Eqs. (22) and fd ae i 
(23) yields Ten = 41g — (2/k)Ro — 2h; (31) 


. For purposes of computation it is more convenient 

-2 dx (22) +4 break the vectors into their real and imaginary 
parts. 


Aw = 2[2 — (j/k)|So — 2S; (24) 4+ QR, +i ~R, (32) 
_%\5 _ 
= + +2/ si la = (3) 
Au = Fs, - 2) Si+S: (26) Ry = lo (34) 
Ay, = 22 — J 2(1 j 
1 2 1 
=2(5 +5) Rot Rt 
* There is little to choose between the mid-chord and the lead- 2 2 1 
ing edge as far as simplification is concerned, and both are su- Pee - i I; +a Rs (36) 


Perior to the quarter-chord for this purpose. 
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2 1 2 1 


2 2 1 
+ + (37) 


In Appendices A and B the foregoing results are 
given for asymptotically small k and large M, re- 
spectively. 


COMPARISON WITH FOREIGN WORK 


It has been brought to the author’s attention that the | 


present problem has been solved by Borbely® in Ger- 
many, Possio’ in Italy, and by the British in works as 
yet unclassified. 

Borbely’s approach is essentially an operational 
solution (via contour integrals) of the equations. His 
results are stated for the lift and diving and moment 
about the leading edge due to plunging and pitching 
about an arbitrary point (on the x axis) and agree 
with the foregoing when reduced to a common pitching 
axis. Since the Green’s function [Eqs. (12) and (16)] 
in the present analysis is stated for an arbitrary velocity 
distribution, it should be more useful than Borbely’s 
results for extension to motions other than plunging 
and pitching. No numerical results were published 
by Borbely. 

Possio’ developed the Green’s function in a rather 
cumbersome series in 6? which appears to be less suited 
to numerical work [even though the J, of Eq. (28) 
may best be evaluated in series form for small values 
of 6 = 2xM]. Possio gives some numerical results 
for M = 1.2, 1.4, 1.6, and 2 and for k = 0 to 0.25. 


NUMERICAL WORK 


Probably the most direct way to evaluate Eq. (28) 
is by graphical. (or numerical) integration, but the 
integrand oscillates so violently for large x and n = 
1, 2, 3 that such a method is too inaccurate. A second 
method is to utilize a power series development of 


Jo(v); since* 
= 2m = sat 
Joo) = an = 
Tp = (40) 


Integration by parts yields 
Leo 


To = = 2xM 
M. 


(41) 


Another approach is effected by introducing the 
expansion‘ 


1947 


1 


m 
Mu ( M ) Mu\” 
m=0 
in Eq. (28), after which repeated partial integration 
yields 


So(8, k) = [Jm(0) + 
GI m+1(8)] = Ro + jlo (43) 


The latter result appears to converge more rapidly 
than Eq. (39) for large 6, although it is useful only for 
integral values of 6.* In order to make effective use 
of Eq. (43), it is necessary to obtain a recursion formula 
for S,; partial integration of Eq. (28) yields y 
M 6 

0 
Borbely gives an analogous result to Eq. (44), although 
the factor 7 of J,-1 was erroneously omitted, at least 
in the R.T.P. translation.® 

For numerical work it is convenient to split Eqs. 
(43) and (44) into real and imaginary parts to give 


— 2n)S,-1 + 2(1 — n)? (44) 


k) = [Jm(8) cos + 


sin 6] (45) 


km 
k) = [—Jm(0) sin + 


Jm+1(8) cos 6} (46) 


= M+ Ih (47) 

=N-—R (48) 

kR, = 2M — P + 31, + (2Ro/6) (49) 
= 2N + Q — 3Ri + (21/8) (50) 
kR; = 4M — 4P + + (8R,/8) (51) 


kI; = 4N + 40 — 5R. + (8h/6) * (52) 


P = (4/mr)Jo(8/M) (cos 6/6) (55) 

Q = (4/r)Jo(6/M)(sin 6/6) (56) 


* At least using Jahnke and Emde;* more extensive tables of 
Bessel functions may now be available. 
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(44) Fic. 4. k*Rea vs. k 
10ugh Fic. 2. k?Rea vs. k. 
(45) / 16 
VA ES 
\ 
4 Z| \ \ 
|| 
(49) 
(50) | 0 4 ar\\ 
| 
0 4 08 k 
| Fic. 5. k*Ica vs. k. 
(52) Fic. 3. k*Zca vs. k. 
(53) to keep them from rising too rapidly and to dampen 
The case of small @ is treated in Appendix A. For otherwise rather violent oscillations for small values 


large @ the asymptotic representation of the Bessel of k. The results agree with those of Possio for 
(54) functions has been tried in Eqs. (28) and (43), but it small k. ‘ 

does not lead to results that are any more useful for In carrying out the computations, the results of 
(55) computation. Appendix A were used for small k, while Eqs. (43) and 
(56) Results have been computed and plotted in Figs. (56) were used for integral values of 6. The con- 

2 through 9 for M = 1.2, ~/2 (labeled 1.4), 1.6, and vergence of Eq. (43) was rather slow (about six terms 
= 2.0 and for k = 0 to 1.5. It was found more con- required) for small M and large k, which factor led 

venient to plot the coefficients multiplied by k? in order to the choice of the upper limit of k = 1.5. 7 
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ai \2 k 


Nn 


Fic. 8. k*Raa vs. k. ’ 
Fic. 6. k?Ray vs. k. Ry 

f | 72 2 2 
Ih = —0/x I! 

AN n= as) 


(Al) 


Mal2 Q, \ 8f1 62/1 1 


where 6° is retained due to the 1/k? factors in Eqs. ly’ 


ig 08 7 Z. (29) through (37). 
Substituting Eq. (Al) into Eqs. (29)-(37) yields 
‘the asymptotic limits. 


Fic. 7. vs. k. 


Rea © — 1) 
CoMPARISON WITH SuUBSONIC RESULTS In = —4/ak 
S In order to facilitate comparison of the foregoing Rea = 4/ak* 
: results with their subsonic counterparts, the incom- Tew = —4/x(M? — Ik 
q pressible results (M = 0) of Theodorsen! and the (A2) : 
Re = 4/3e(M? — 1) I, 
compressible results of Possio* (for M = 0.7) have 
been plotted wherever feasible in Figs. 2 through 9. Tan = —2M*k/3x(M? — 1)? 
For M = 0, Ian = 0. Ria = 8M?(2M? + 1)/34(M? — 
Ten = 4(2 — M?)/3e(M? — 


APPENDIX A.—SMALL VALUES OF k 


The limiting values of R, and S, for small 6 are ~ For more general motions Eq. (12) may be asymp- The 
obtained from Eq. (28) as totically expressed as be 


» 

le 

7 
/\ x \ R= + — 
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ate, ) = 1 - @ - x 
(3+ Ga) ak Hee 
7, 
APPENDIX B.—LARGE cF M / 
For large values of M, Eq. (12) is asymptotic to /; V / 2 
while S, becomes YW 
Jo M? 
(B2) k 
TOn = On JRIn+1 nt2 a8 4 
Evaluating Eq. (B3) yields a2 
R! = s/k ) Fic. 9. vs. k. 
Ri! = —(1/k*) + (c/k*) + (2s/k) mR, = ~ ie Ry+2’ — | 
Ry’ = (4c/k*) + (2/k)[2 — (1/k*)]s (BS) 
2 2 1 4s = I,’ - + | = 
1 
(Al) I,’ =I," 
6 6 1 4 3 M* 
Substituting Eqs. (B4) into Eqs. (B5) yields 
ke) In!’ = —(1/2k) + [k + (1/2k)]c + 
24 = [2k + (1/k)]c + [1 — (1/2k?)]s 
elds 16 (2 3 ) I,’’ = 4ke 
Re) 6 > (B6) 
is (2 2 3 
80 ) 
I) = 2k) c+ 6 3 
ks 8 Substituting Eqs. (B4)-(B6) into Eqs. (29)-(37) yields 
¢ = cos 2k,s = sin 2k = —(1/M*)(2s/k) (B7) 4 
The R and J are calculated from Eqs. (B2) which may 
be rewritten 
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Coupled Vibrations of an Elastically Mounted 
Rigid Body with One Plane of Symmetry 


R. J. DUNHOLTERt 


Unwersity of Cincinnati 


ABSTRACT 


An elastically mounted in-line aircraft engine is treated as a 
rigid body with a vertical plane of symmetry. The symmetry 
permits the six equations of motion to be resolved into two inde- 
pendent groups of three each. By using dimensionally homogene- 
ous variables and the methods of vector-matrix analysis, the three 
scalar equations in each group are reduced to one. All the inertial 
and stiffness properties of the engine and its mounting system are 
contained in a single matrix coefficient called the frequency 
matrix. The static displacement, natural frequencies, and, the 
dynamic magnification factors are expressed as operations on the 
frequency matrix. A geometric interpretation is given for the 
modes of motion and the relative amplitude vectors. 


INTRODUCTION 


A* IN-LINE AIRCRAFT ENGINE can be assumed to 
have a vertical plane of symmetry. It is custo- 
mary to mount the engine at four points lying in a plane 
normal to the plane of symmetry and containing the 
crankshaft centerline. In general, the center of mass of 
the engine will not lie in the plane of the mounting 
points, nor will the longitudinal principal central axis be 
parallel to this plane. Thus, there will be coupling 
terms in all six of the equations of motion. By a suitable 
choice of the elastic properties of the mounting devices 
it is possible to make some of the coupling terms vanish. 
It is not practicable to make them all vanish. It is 
desirable to have a general solution describing the modes 
of motion and giving the natural frequencies, the dis- 
placements, and the dynamic magnification factors. 
This problem has been discussed by von Schlippe.' 
The solution presented here is based on a generalization 
of the simple linear vibration equation 


¥+ w,’x = 0 (a) 


and is suitable for routine calculation. 

Only the case of a propeller at rest will be considered 
here. Maximum simplicity ot the equations of motion 
is achieved by using the principal central axes for refer- 
ence. The positive directions of the three linear and 
three angular coordinates are shown in Fig. 1. 

Because of the plane of symmetry (XZ) the equa- 
tions of motion consist of two independent sets of three 
each. The equations containing the coordinates x, 0, z 
describe the symmetric modes of motion, and the equa- 
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tions containing the coordinates ¢, y, y describe the 
antisymmetric modes of motion. 


SYMMETRIC MODES 


Undamped Free Vibrations 
The scalar equations of motion are 
mi + + Co + C2 = 0 
Tp + Corx + Coo + Coz = 0 (1) 
me + Cx + Co + Cyt = 
where m is the total mass, J, is the mass moment of 
inertia with respect to the Y axis, and C;,, C,,, etc., 
are the stiffness coefficients for the complete engine 
mount. The set of Eq. (1) can be written in the vector- 
matrix form 


=0 (2) 

where the displacement vector r’ is 
r’ = [x, 6, 2] (3) 
the inertia matrix J is 


J, 0 (4) 
00m 
and the stiffness Matrix V is 
Ces Cu Cus 
(Cor Co Cr (5) 
Cuz Cu Cus 
The displacement vector 1’ is not dimensionally 
homogeneous. The vector 


r= = = 


[Vmx, Vmz] (6) 


is dimensionally homogeneous and will be called the 
d.h. displacement vector. Let Eq. (2) be multiplied 
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by I-'’? and the unit matrix U = I~’/- I’ inserted as 
shown in Eq. (7). 


4+ = 0 (7) 
By calculation, we find that 


m m 
C, C, 

Oz “00 62 

m Vml, m 


The matrix of Eq. (8) will be called the natural fre- 
quency matrix Q and the elements designated as shown 


in Eq. (9). 


Q = | wz” woe? 9,” (9) 
@22” 


The elements in the principal diagonal of the matrix 2 
are necessarily positive, but the other elements may be 
negative. 

Eq. (7) can now be written in the form 


7T+Q-r=0 (10) 


which is a generalization of Eq. (a). A particular solu- 
tion of Eq. (10) has the form 


= 7% sin wt 


(11) 


where 7% is the d.h. amplitude vector. Substituting r 


of Eq. (11) in Eq. (10) one obtains 


= (12) 


The natural frequency matrix 2 is symmetric. It is 
well known from matrix theory that Eq. (12) will then 
have, in general, three distinct solutions. That is, 
there are three sets of w? and 1% which satisfy Eq. (12). 
The directions of the solutions 7 are unique but the 
magnitudes are not. 


Fic. 2. The dimensionally homogeneous amplitude vector 7. 
Symmetric modes. 
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The values of w’ are the roots of the natural frequency 
equation 


(w?)? — + — 2 = 0 (13) 


where 2, is the sum of the elements in the principal dj- 
agonal of Q, {2 is the sum of the cofactors in the principal 
diagonal of Q, and Q; is the determinant of Q. 

With each w,?, i = 1, 2, 3, there is associated a value of 
(7%), 4 = 1, 2, 3, given by the formula 


(To) 4 = (Q i= 1, 2,3 (14) 


where v is an arbitrary vector and (Q — 00°U) aay i is the 
adjoint matrix of (Q — w,?U). 

Taking v = j = [0, 1, 0], the components of (1), 
are found to be 


V = + — 

+ wat} (15) 
m( 20): = + — 


Furthermore, it is known that the solutions for (1%), 
of Eq. (12) are mutually perpendicular. Thus, 


(To) 4° (To); = 0, #=j=1,2,3 (16) 


Eq. (16) is the orthogonality condition on the natural 
mode amplitudes. Fig. 2 shows a geometric interpreta- 
tion of the natural mode d.h. amplitude vector 1%. 


The projection of 7 on the 6 axis-(Y), that is OR, 
gives the direction of the axis of rotation. The magni- 


tude of OR is +/J, times the amplitude of the angular 
displacément %. The projection of 7 on the XZ plane, 


ij, 


that is OD, gives the direction of the linear displacement 


of the center of mass. The components of OD on the 
X and Z axes are ~/m times the amplitude of the x 
displacement x and +/m times the amplitude of the z 
displacement 2, respectively. 


It is apparent from Fig. 2 that the motion of the 
body is a rotation about the Y axis followed by a trans- 
lation of the center of mass in the XZ plane. This 
motion is equivalent to a pure rotation of the body 
about an axis normal to the XZ plane whose trace has 
the coordinates 


(px): 
+=1,2,3 (17) 
— (x0) (80) 


The position vector p; of the trace point P; can be 
expressed in the vector form 


y 
i= 

mj 
where j is a unit vector in the direction of the Y axis. 
If the trace points of the axes of rotation are connected 
to form a triangle as shown in Fig. 3, then the center of 


mass lies at the point of intersection of the altitudes. 
It is easily shown, for example, that 


(2) (An, 


i=1,2,3 (18) 
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P1‘(p2 — ps) = 0 
by using Eq. (18) and Eq. (16). Thus the above state- 
ment is proved. This geometric property can be used 
as a check on the numerical calculations. It follows 
from Eq. (18) that, since the directions of (1%); are 
unique, the position vectors p; are unique. 


Static Displacement 
It is apparent from Eq. (2) that the generalized vec- 
tor force P’ acting on the body due to a displacement 
vector r’ is 
P’ = -v-r’ (19) 
If F’ is an applied force, then F’ = —P’ for equilibrium 
and 
F’ = W-r’ (20) 
Eq. (20) can be expressed in a dimensionally homogene- 
ous form involving the natural frequency matrix Q as 
follows: 


17. 
Thus, 
F=aQ-r (21) 


where 


1 1 1 
Vm "Vl, Vm 
is a d.h. force vector and r is the d.h. displacement vec- 


tor. 
Solving Eq. (21) for r, one obtains 


r= Q-'-F (22) 


Thus, the reciprocal of the natural frequency matrix is 
the dimensionally homogeneous displacement matrix. 
Eq. (22) will be used to obtain the dynamic magnifica- 
tion factors. The displacements alone are more easily 
calculated from Eq. (20), which gives 


= (23) 
Forced Vibrations 
A generalized force of the form 
F’ = F,’ sin wi (24) 
where 


Fy’ = [F:, M,, F,] 


will excite only the symmetric modes. The equation of 
motion is 


+ = F’ (25) 
The dimensionally homogeneous form is 
r+Q-r=F (26) 


Substituting 


/ 


/ 


Fic. 3. The traces P; of the axes of rotation. Symmetric modes. 


~ 


r = 7% sin wt 
in Eq. (26), one obtains - 
+ Q-% = Fo 
(Q — wU)-m = Fy (27) 
where U is the unit matrix. 
The solution of Eq. (27) is 

TOaynamie = (Q — w?U)—'- Fy (28) 

The factor Q-Q-! = Vis inserted in Eq. (28) as follows: 
= (Q — F, 

(Q — 


= 


or 


Thus, 


where M = (U — w’Q-')—' is the dynamic magnifica- 
tion matrix for the d.h. displacements. Let the ele- 
ments of the matrix M be designated M,,, Mz, Mzz, 
etc. The existence of an element M,, indicates that 
there will be a dynamic d.h. displacement in the x 
direction (linear) corresponding to a static d.h. dis- 
placement in the @ direction (angular displacement 
about Y). The ratio of these displacements is the 
number MM. 20+ 

The magnification matrix can be interpreted also in 
another way. Since the natural mode d.h. amplitude 
vectors (7%); 7 = 1, 2, 3, are mutually perpendicular, 
they can be used as reference axes. Let uj, Us, U2 be 
unit vectors in the directions of (1)1, (To)2, (To)s, respec- 
tively. Let these normal axes be used in Eq. (26). 
The components of the d.h. force amplitude vector Fo 
are components in the directions u;, Uz, Us. Since 
U, Uz, Us represent the three natural modes of motion 
it follows that, if Fp has a component on u,, for example, 
then that particular mode of motion will be excited. 
Corresponding to the component of Fp on u; there is a 
static deflection of the body which, in general, involves 
all three coordinates x, 8,2. Depending upon the excita- 
tion frequency , this static deflection in the first nat- 
ural mode will be magnified in the ratio given by the 
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matrix M. For the base vectors u;, Us, U3, the matrix M 
has the form : 


1 
i-wie) ° 
1 
M= 0 0 (30) 
1 
? 


Thus the magnification for the first natural mode is 
1/[1 — (w?/«’)]. 
Let Eq. (29) be multiplied by Q. Thus, 

Q-Toaynamic = Q-M- Tostatic (31) 
Since Q and M are symmetric, one can write 

Noaynamic = Tostatic 
or, using Eq. (21), 

Foaynamic = M- Fostatic (32) 


Thus M is the dynamic magnification matrix for both 
d.h. displacement and d.h. force. 


ANTISYMMETRIC MODES 


The equations for the antisymmetric modes are for- 
mally identical with those for the symmetric modes. 
The d.h. displacement vector and the d.h. amplitude 
force vector that excites these modes are, of course, 
different. 


Undamped Free Vibrations 
The vector-matrix equation is ° 
7+Q-r=0 (33) 
where the d.h. displacement vector is 


r= [ Vie, V my, (34) 
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Fic. 4, The dimensionally homogeneous amplitude vector fo. 
Antisymmetric modes. 
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xX 


Fic. 5. Axis of pure rotation. Antisymmetric modes. 


and the natural frequency matrix is 


Cog Coy Coy 


I, V mi, 
(35) 


(36) 


The components of the natural mode vectors are 
= Wyy + — wo, *wyy? 
= (02)? — + + 
— 


The orthogonality condition on the natural modes of 
motion in expanded form is 


+ men) + = 0,4 (38) 


Fig. 4 shows a geometrical interpretation of a natural 
mode vector. The projection of ™ on the ¢y(XZ) 
plane is not, in general, the true axis of rotation. The 
exception occurs when J, = J,. If the components of 


OR on ¢ and y are divided by +/J, and ¥/I,, respec- 


tively, one obtains the components of OR’, which is 
the true axis of rotation. Thus, in a natural mode of 
motion the body rotates about the axis OR’ through 
the angle +~/(¢o)? + (Yo)? as the center of mass moves 
in the direction of the Y axis a distance yp. The nat- 
ural mode motion can also be described as a pure rota- 
tion about an axis lying in the XZ plane and parallel 


(37) 


to OR’. The position of this axis of pure rotation can 
be determined by the intercepts p, and p, on the X 
and Z axes, respectively (Fig. 5). 
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The total angular rotation is 
dot + Yok 
Then the displacement of the center of mass is 
(Goi + Yok) X (—p.k) = 
(Goi + Yok) X (—pai) = 


From the above equations the intercepts are 


(39) 


The three axes of pure rotation form a triangle 
having the following two properties: (1) Each vertex is 
the position of the resultant inertia force corresponding 
to the rotation of the body about the side opposite; 
and (2) the center of mass lies within the triangle. Fig. 
6 shows the directions of the inertia moments and the 
inertia force for rotation about the given axis. 

It is apparent from Fig. 6 that the coordinates of the 
point C, the position of the resultant inertia force, are 


= 2¥0/MYo —I,/mp, 
(40) 


The equation of the natural mode axis (Fig. 5) is 


2 = + = 1, 2,3 (41) 

If C, lies at the intersection of the axes i, 7 ~ j, then 
T, _ (bite 1 


Clearing Eq. (42) of fractions and substituting the 
values of p, and p, from Eq. (39), one obtains 


This is identical with the orthogonality condition, Eq. 
(38), and thus the supposition is proved. It is apparent 
from Fig. (6) that C and the center of mass CG always 
lie on the same side of the axis of rotation. Since C is a 


+ t AZ (42) 
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Fic. 6. Inertia force and inertia moments. Antisymmetric 


modes. 


vertex of the triangle formed by the axes of rotation, it 
follows that the center of mass CG lies within the tri- 
angle. The line joining a vertex with CG will not, in 
general, be normal to the side opposite. The exception 
occurs when J, = J,. 


Forced Vibration 
A generalized force of the form 
F’ = sin wt 
where 
[M,, Fy, M,) (43) 
will excite only the antisymmetric modes. The d.h. 
force amplitude vector is 
REFERENCE 
von Schlippe, B., Schwingungsberechnung von raumlichen 
Masch ten, Jahrbuch 1937 der Deutschen Luftfahrt- 


forschung, s. II 103. 
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Determination of Airplane Thermometer 
Recovery Factors in Flight 


H. W. SIBERT* 


University of Colorado 


SUMMARY 


A formula, which is apparently new, is derived for the airplane 
thermometer recovery factor as a function of the free-stream 
Mach Number. A semigraphical method, based on this formula, 
is then outlined for determining the thermometer recovery factor 
by flying the airplane at various free-stream Mach Numbers at a 
constant altitude. 


INTRODUCTION 


A THERMOMETER placed in a stream of air moving 
relatively to it will register a higher temperature 
than that of the surrounding air. This is because a 
portion of the air in the vicinity of the thermometer has 
been slowed down relative to the thermometer. In a 
typical airplane thermometer installation for determin- 
ing the free-stream air temperature, the thermometer is 
surrounded by a cylinder that is open at its forward end 
and has several holes in its rear portion. This arrange- 
ment ensures a continuous, but relatively slow, flow 
past the thermometer. Formulas’for obtaining the 
recovery factor for such a thermometer will be derived 
below for a thermometer installation that is entirely 
outside of any boundary layer. 


DERIVATION OF FORMULAS 


Let the thermometer installation be outside of any 
boundary layer, and let 


T, = free-stream absolute temperature 

v, = free-stream speed relative to the thermom- 
eter 

absolute temperature registered by the ther- 
mometer 


vm = speed of air relative to the thermometer and’ 


corresponding to T 
k =thermometer recovery factor = 1 — 
If there is no shock wave in front of the thermometer, 
the flow will be isentropic and 


+ CT, + C,T 


in which C, is the specific heat at constant pressure. 
With a shock wave in front of the thermometer, the 
flow will still be isentropic on either side of the shock 
wave. Hence, 
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+ Gz, = C,T a 
+ = CT 2 
where and T are, respectively, the stagnation tem- 
peratures on the upstream and downstream sides of the 
shock wave. It is well known that 7, = T,. Hence, 
the equation 
+ CT, = + GT (1) 
is valid with or without a shock wave in front of the 
thermometer. From this equation 
— T,) = 1/2,? = (2) 
In terms of the ratio of specific heats y and the gas 
constant R, 
C, = yR/(y — 1) (3) 
Moreover, the speed of sound a is defined by the rela- 
tion 


a? = yRT (4) 
From these last three equations : 
(Ta — T,)/T, = — 1)ko,?/a,? (5) 


in which a, is the free-stream speed of sound. For y = 
1.4, Eq. (5) becomes 
Ty, = T,(1 + 0.2kM,?) (6) 
in which M, is the free-stream Mach Number. Note 
that 
T =t+ 459.4 (7) 
where ¢ is in degrees Fahrenheit. From the last two 
equations 
ty = t, + 0.2(t, + 459.4)kM,? (8) 


DETERMINATION OF RECOVERY FACTOR IN FLIGHT 


By means of Eq. (8), & can be determined by flying 
the airplane at various values of M, at a constant alti- 
tude (so that ¢, will remain essentially constant). A 
curve of ft, versus M,? is then plotted, and the slope m 
of this curve is determined for various values of M,. 
Since ¢, is essentially constant, it follows from Eq. (8) 
that the value of k at each of these values of M, is given 
by the relations 

0.2(t, + 459.4)k = 


k = 5m/(t, + 459.4) 
(Continued on page 367) 


(9) 
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Effect of Wind Gradient on Glide and Climb 


BERNARD ETKIN* 
Unwersity of Toronto 


SUMMARY 


The problem of an aircraft flying through a head wind that 
increases with altitude is examined. The equations of motion 
are set up and solved for flight at constant air speed in a wind 
that varies linearly. It is found that to fly under these condi- 
tions the angle of attack must be quickly adjusted to a steady 
value and that the flight is subsequently at constant attitude, 
angle of attack, and elevator angle. The path is parabolic, the 
glide being extended and the climb contracted. Graphs are 
given of glide characteristics for a typical light aircraft. 


INTRODUCTION 


[ IS WELL KNOWN that wind velocity generally in- 
creases with altitude at the lower levels. The effect 
of this gradient wind on an aircraft flying through it 
does not seem to have been investigated, although on 
one occasion incorrect conclusions have been published.' 

An approximation to the steady-state variation of 
wind speed with height is given by the Ekman spiral 
(reference 2, Chapt. X). For typical cases the graph 
does not deviate much from a straight line up to about 
1,000 ft. 

Experimental data also show that the variation near 
the ground is often nearly linear and that the wind at 
1,000 ft. may be as much as ten times that at the ground 
(reference 3, pp. 256, 244). In the present paper a 
solution is given to the motion of an aircraft climbing 
or gliding in a linearly varying wind. 


EQUATIONS OF MOTION 


The equations are written and solved for the gliding 
case first, the analysis for climbing following directly. 
Referring to Fig. 1, the equations of motion are 


gS (C, sin — Cp cos g) = 
gS (C, cos ¢ + Cp sin ¢) -—- W= my 
gScCy = —I6 


where 


= '/)V? 

= air density 

= wing area 

= mean chord 

= weight of aircraft 

coefficients of lift, drag, and moment 
=, pitching moment of inertia 

= mass of aircraft 

= d*x/dt* 

d*y/dt? 


Se: & 3 HrR 
ll 
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The component velocities and accelerations are given 
by 


y = —Vsing 
y = —Vsing — 
x = Vcecose—w 


x= V cos ¢ — V sin g¢ — wy = 
V cos g — V sing (¢ — w,) 
where 


w = W+ wy = head wind 
Wo head wind at ground 
Wy dw/dy 


Thus the equations of motion become 


V cos g — Vsin — w,) = 


2 (C, sing — Cpcos¢) (1) 


S 
V sin g + V cos gg = g — * (C, cos ¢ + Cp sin ¢) 


: (2) 
I6 = —CyqSc (3) 


Of the quantities that occur in these equations, C, 
and Cp are functions of a; Cy isa function (when damp- 
ing is considered) of a, 6, 6, where 4 is the elevator angle; 
and ¢ isa function of a,@. Thus the equations contain 
the four independent unknowns I, a, 0, and 6. One of 
these may be assigned an arbitrary value, either con- 
stant or variable with time, and a solution may be 
sought for the remaining three. In general, the equa- 
tions are nonlinear and not soluble by ordinary means. 
Fortunately, however, the case of most interest is that 
of flight at constant air speed, and this case can be 
solved explicitly under the assumption of small angles, 
which is a good approximation. 5 


HORIZONTAL 


ZERO LIFT DIR'N 
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FIXED AXES 
re) x 
Fic. 1. V = airspeed vector. w = wind vector. 
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ASYMPTOTE 
10 
| Rez 


TIME - SECS. 


Fic. 2. Example of variation of a with time for constant air 
speed. 


SOLUTION FOR CONSTANT AIR SPEED 


With V = constant and sin g = ¢, cos g = 1, the 
equations become 

—Vo(e — wy) = (gS/m)(Cre — Cy) (4) 

Ve = g — (gS/m)(C, + Cog) (5) 

= —CyqgSc (6) 

Since Eqs. (4) and (5) contain only the implicit un- 

knowns a and 6, they may be solved without regard for 


Eq. (6). The solutions for a and @ used in Eq. (6) will 


then serve to determine the elevator angle 6. 
Eliminating ¢ from Eqs. (4) and (5) yields, on ne- 
glecting ¢’, 
¢ = gSCp/m(g — w,V) 

or 

g=AC (7) 
where A = constant. Substituting Eq. (7) in Eq. 
(4) gives 

—VA¢g + VAw, = (gS/m)(AC, — 1) (8) 


Writing C, = aa and Cp = Coo + Ka’, where a, K, 
Cpo are constants, and noting that 


= A(dC)/dt) = 2KAa& 
Eq. (8) becomes 
pVS | pV Sa 
aL2KA 


where 


P 
4G 
4mK A? 4mKA a 2 ©) 


P = (w,/2KA) + (pVS/4mKA?) 
Q = pVSa/4mKA 


The solution to Eq. (9) is given by 


a 


or 


a P — Qao 

i= 2 

+ (P/Q) loge Oa 

It is noted from Eq. (9) that when a = P/Q, then 

& = 0 and a steady state exists. When the angle of at- 
tack differs initially from this value, Eq. (10) gives the 


(10) 


variation with time which must exist if the air speed is 


to remain constant. As illustrated in Fig. 2, the angle 
must be brought very nearly to the asymptotic value 
a’ = P/Q within a short time, of the order of 1 sec. 
The computation of Fig. 2 is for a light airplane gliding 
at 60 m.p.h. through a wind gradient w, = 0.044 ft. 
per sec. per ft. The airplane has the following charac- 
teristics: 


W = 1,200 lbs. 

S = 170 sqft. 

Cp = 0.030 + 1.08a? 
CL = 4.30a 


The critical value of a is a’ = 10.20°. For an arbi- 
trary initial value a9 = 7.96°, the variation is as shown. 
Since the angle of attack must so rapidly be adjusted 
to the value a’, it may be considered that the constant 
air-speed glide is also at constant a. Eq. (7) yields that 
yg must also be constant; hence, 6 = constant. From 
Eq. (6) it follows that C,, and 6 are constant as well— 
that is, no movement of the stick is required to main- 
tain the steady glide at constant air speed. 


GLIDE PATH 


In the steady-state condition the velocity and posi- 
tion of the airplane are given by 


p= -Vsing 
y=m — Vtsing (11) 


x = Veosy — w= Vecosg — (w + wy) = V cos 
¢ — [wo + — Ve sin 


therefore 


x = Vt cos — (Ww + + sin = 
(V cos ¢ — w)t + 1/swyV sin gl? (12) 


where w, = head wind at initial height 4. 

The path is given by Eqs. (11) and (12) and is para- 
bolic, as though the aircraft were a projectile in a hori- 
zontal gravity field of acceleration w,V sin ¢. The 
path is drawn in Fig. 3 for the same aircraft chosen for 
illustration above. It will be noted that the devi- 
ation from the initial tangent is considerable in this 
case. 


“OvERSHOOT’’ DISTANCE 


The airplane overshoots the initial tangent by the 
distance s = '/.w,V#? sin g. This quantity has been 
computed for a 1,000-ft. descent at three wind gradients. 
The results, given in Fig. 4, are for the same aircraft 
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FLIGHT PATH jor 


Gat 
HORIZ. DIST. +1000 FT. MPH 


Fic. 3. Glide path of typical - airplane. Air speed = 60 
m.p.h. 
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EFFECT OF WIND GRADIENT ON GLIDE AND CLIMB 


he angle as used above. It will be observed that the wind gradi- ‘ARIES ECE ed 
Ic value ent may be responsible for considerable overshoot, / INO GRADIENT 
bf 1 sec. perhaps more than the pilot can compensate for, 
P gliding and provides an additional hazard in forced land- 8 Bok 
044 ft. ings. 
charac- 
EFFECT ON CLIMB 
Ra: 
The problem of the climbing aircraft differs from the HA 
glide in that the angles ¢ and @ are negative and that a ee 
there is an additional force, the thrust 7. The addi- 5 1 2 ——— es 
tion of a constant thrust to a steady-state solution will re} as ot Pee 
n arbi- yield another steady state. Hence, the climb at con- 2 ane 
shown. stant air speed will also be at constant a, 0, g. The wi 
A ath will be a parabola as before, with the convexit yA 
in the same sense, so that the horizontal distance to 
AIRSPEED - MPH 
ds that reach a given height is contracted as compared with the ; 
From case of no wind gradient. Fic. 4. Overshoot distance vs. air speed and wind gradient for a 
1 typical light airplane in 1,000-ft. descent. 
well— 
, main EFFECT OF AIRPLANE CHARACTERISTICS 
The distortion of the paths in climb and glide is esiana ict 
greatest when the rate of descent is small and when the ‘ Lloyd, Sloan, Best Speed for Climb, Can. Aviation, October, 
wind velocity is large relative to the airplane speed. 198- | 
i ‘ 4 5 ? Haurwitz, Dynamic Meteorology; McGraw-Hill Book Com- 
1 posi- Thus, the aircraft most affected would be machines with cnnty, ths... Haul anaes 
low wing loading and flat glides. Fast machines with ~ : Byers, Synoptic and Aeronautical Meteorology; McGraw- 
steep glides would be less disturbed. Hill Book Company, Inc., New York, 1937. 
(11) 
V cos 
Determination of Airplane Thermometer Recovery Factors in Flight 
(Continued from page 364) 
= 
(12) The value of ¢, to use in Eq. (9) can be determined 
directly by means of a thermometer attached to a 
sounding balloon, or it can be found as the /-intercept 
para g 
hori- obtained by extending the ¢,, versus M, curve to M, = 
The 0 [see Eq. (8)]. 
n for It should be noted that if temperatures are measured 
devi- in degrees Centigrade, Eqs. (8) and (9) still hold pro- 
- this vided the number 459.4 is changed to 273. 
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